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I.  INTRODUCTION 


Rate  is  a  concept  common  to  almost  all  scientific  and  technical  disciplines.  The  Oxford 
Desk  Dictionary  [Abate  1997]  defines  rate  as  a  “numerical  proportion  expressed  in  units” 
and  as  a  “pace  of  movement  or  change.”  These  definitions  are  not  very  helpful,  but  they  do 
convey  the  deceptive  simplicity  of  the  idea. 

Despite  this  almost  pervasiveness,  there  is  little  in  the  way  of  a  common  theory  of  rates. 
In  particular,  the  various  formulations  of  rates  are  often  imprecise  and  ad  hoc.  Further,  the 
rate  differential  equation,  an  inseparable  part  of  the  basic  concept,  is  almost  always  a  deus  ex 
machina.  As  a  result,  rate  formulations  are  often  arcane  and  less  accurate  than  they  could 
be. 


Our  purpose  here  is  to  investigate  and  elaborate  the  nucleus  of  just  such  a  general  theory 
of  rates.  This  theory  is  new  and  incomplete,  but  it  offers  great  potential  for  the  understanding 
and  application  of  rates. 

This  body  of  methodology  was  first  introduced  in  [Fowler  2001]  and  then  elaborated  in 
[Fowler  2004]  to  expand  the  nature  of  attrition  rate  differential  equations. 

II.  WHAT’S  A  RATE? 

Regardless  of  whether  we  are  describing  the  collisions  of  photons  with  other  particles,  the 
chemical  interaction  of  molecules,  or  the  populations  of  sexually  or  asexually  reproducing 
animals,  we  are  dealing  with  events  that  are  well  localized  in  time;  that  is,  they  take  less 
time  to  occur  than  the  period  of  observation.  To  some  degree,  we  may  view  these  events  as 
occurring  impulsively  or  instantaneously,  although  we  shall  almost  immediately  generalize 
the  theory  to  accommodate  events  that  have  a  finite  duration. 

This  leads  us  to  a  third  definition  of  rate:  rate  is  a  means  of  continuously  representing 
the  change  in  quantity  of  a  population  of  discrete  events  or  other  observables.  By  this  we 
mean  that  a  rate  is  a  mapping  of  discrete  events  into  a  continuous  representation. 

Actually,  we  want  to  also  consider  the  opposite  as  well,  the  mapping  of  a  continuous 
process  into  discreteness.  For  simplicity  at  this  point,  we  shall  assume  this  complementary 
aspect  is  merely  a  matter  of  piecewise  integration,  however  deceptively. 

Simply  put,  however,  a  rate  is  a  mathematical  expression  of  how  many  things  take  place 
in  a  unit  of  time.  As  such,  it  is  a  change  of  representation  from  a  collective  discrete  rep¬ 
resentation  of  when  each  of  the  things  occur  to  a  continuous  representation.  Ideally  this 
transformation  reproduces  the  occurrence  times  of  the  events  accurately,  although  the  very 
concept  of  a  rate  implies  some  sort  of  smoothing  or  averaging  of  event  times. 

We  advance  a  set  of  propositions  to  lay  the  groundwork  for  a  general  theory  of  rates: 
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Proposition  1.  A  rate  is  a  mapping  of  a  collection  or  series  of  connected  but  discrete  events 
or  observations  into  a  continuous  counting  representation. 


This  is  essentially  a  restatement  of  the  basic  definition  above  in  somewhat  more  math¬ 
ematical  terms.  The  key  point  is  that  we  are  starting  with  discrete,  that  is,  punctuated 
in  position  (normally  time),  events  that  we  transform  into  some  continuous  representation. 
The  continuous  representation  is  an  enumeration  of  the  events  even  though  it  is  no  longer 
integer. 


Proposition  2.  Observable  events  consist  of  at  least  two  measurable  quantities:  the  imme¬ 
diate  intensity  arising  from  the  occurrence  of  the  event  (which  may  be  just  the  occurrence) 
and  the  position  of  the  event. 


Each  of  these  punctuated  events  has  a  magnitude  and  a  position.  That  is,  we  know  how 
big  some  aspect  of  the  event  is.  Assuming  the  position  is  temporal,  we  know  when  the  event 
occurs. 

At  this  point,  we  have  not  yet  addressed  whether  the  duration  of  the  event  is  finite  (as 
opposed  to  infinitesimal.)  We  further  assume  that  these  two  quantities  are  either  completely 
independent  of  each  other,  or  if  they  are  dependent,  we  have  a  model  of  the  nature  of  the 
dependency. 


Proposition  3.  Events  occur  as  the  effect  of  some  cause. 


This  is  not  some  supernatural  phenomenon.  We  assume  that  the  basic  scientific  relation 
of  cause  and  effect  appertain.  We  further  assume  that  we  know  this  relationship. 


Proposition  4.  Events  may  be  caused  by  single  or  multiple  agents. 


We  assume  that  the  causality  can  be  systematized;  that  it  is  not  holistic.  Specifically,  we 
assume  that  the  cause  agents  may  be  aggregable. 


Proposition  5.  Events  may  be  described  by  some  mathematical  formalism. 


The  mechanics  of  the  events  has  some  mathematical  description.  This  follows  largely 
from  the  proposition  of  cause  and  effect,  but  it  also  says  that  things  are  fairly  well  behaved. 


Proposition  6.  The  mechanism  that  causes  events  may  be  singular  or  repetitive. 
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This  simply  says  that  cause  agents  may  cause  one  event  or  they  may  cause  several  in 
succession. 


Proposition  7.  Cause  agents  are  independent  and  may  be  identical,  but  their  association 
in  a  rate  process  has  a  collective  effect. 


This  lets  us  treat  the  agents  as  if  they  are  identical  but  independent.  It  also  acknowledges 
that  they  have  a  different  effect  when  they  band  together.  Recall  Napoleon’s  dictum  ’’There 
is  a  quality  in  quantity.” 


Proposition  8.  Events  occur  in  a  discrete  manner  although  not  necessarily  limited  to  an 
infinitesimal  change  in  position. 


Although  the  events  are  punctuated,  they  may  not  occur  instantaneously. 

Our  intent  is  captured  in  Figure  1  where  we  depict  two  discrete  event  trajectories  (they 
look  rather  like  stairs,)  and  two  equivalent  continuous  representations.  In  this  case,  these 
are  taken  from  the  application  of  General  Rate  Theory  to  attrition. 


Figure  1.  Discrete  Event  and  Continuous  Representation  Trajectories 


Our  approach  is  to  first  briefly  present  the  pieces  of  mathematics  needed  to  develop 
rate  theory.  We  then  use  these  pieces  to  first  develop  the  continuous  representations  of  the 
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discrete  events  and  then  the  general  rate  differential  equations.  Finally,  we  present  several 
examples  of  the  application  of  the  theory. 

III.  FUNCTIONAL  THEORY 


Our  starting  point  in  developing  this  General  Theory  of  Rates  is  what  is  often  called 
Functional  Theory  [Stakgold  1998],  although  it  is  sometimes  called  the  Theory  of  Generalized 
Functions  [Erdelyi  1962]  or  Distribution  Theory.  One  of  the  most  common  uses  of  this  theory 
is  the  development  of  Green  Functions.  A  Green  Function  is  a  representation  of  the  behavior 
of  a  mechanical  system  under  the  imposition  of  an  impulsive  “force."  Our  approach  is  similar 
in  that  we  are  dealing  with  impulsive  events.  The  difference  is  that  the  Green  Function,  via 
the  vehicle  of  the  properties  of  the  Dirac  delta  function,  is  used  to  describe  the  behavior  of  a 
mechanical  system  under  the  imposition  of  a  general  “force, "  while  our  interest  in  a  General 
Theory  of  Rates  is  finding  one  or  more  means  of  transforming  a  general  collection  of  Dirac 
delta  functions  into  a  continuous  representation. 


The  Dirac  or  strong  delta  function,  5  (t  —  t').  is  not  a  function  in  the  normal  sense.  In 
particular,  it  has  the  behavior 


(  0,  t  <  t' 

S  (t  —  t')  —  <  1,  t  =  t' 

{  0 .t>? 

One  then  naturally  asks  the  question,  If  the  integral 

tf  ~ \~£ 

J  5  (t  —  t')  dt'  =  1,  e  >  0, 

t'  —  £ 


then  for  what  value  of  s'  <  e  does  the  equivalent  integral 

t'+e' 

J  5  (t  —  t')  dt'  =  / 

t’-e’ 


(1) 


(2) 


(3) 


have  a  value  of  0  <  /  <  1?  This  is  clearly  conceivable,  and  our  intuition  tells  us  that 
such  a  situation  must  exist.  The  methodology  or  formalism  for  describing  this  is  given  by 

Functional  Theory. 


Delta  functions  have  functional  representations  as  Gaussians  (among  other  forms)  in  a 
weak  sense.  By  weak,  we  mean  that  if  we  consider 


g  (x  :  x,  a)  =  _^e-(*-*)2 /2<^ 


V/27T( 


TUT 


which  we  know  has  unit  integral,  then 


lim  g  (x  :  x,  a)  — >  5  {x  —  x) 

cr—>0 


(4) 

(5) 
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in  the  sense  of  Equation  (1). 


We  may  extend  this  idea  even  to  the  simple  Negative  Exponential  Distribution  (NED), 

pit :  a)  —  a  e~at.  (6) 

If  we  take  the  limit  a  — >  oo,  then 

p((:“^°o)={  M>0  (7> 

This  is  again  a  delta  function.  We  shall  refer  to  this  behavior  as  a  weak  delta  function  (or 
functional)  in  contrast  with  a  classical  or  strict  sense. 

The  basic  idea  of  functional  (or  distribution)  theory  is  that  instead  of  describing  a  function 
f  (x)  with  a  value  at  every  point  x,  we  use  the  real  number  f  f  (x)  (f>  (x)  dx  for  every  0 

Rn 

belonging  to  a  class  of  accessory  functions.  In  this  representation,  /  is  a  functional  on  the 
class  of  accessory  functions.  For  developing  rate  theory,  the  class  of  accessory  functions  is 
comprised  of  the  individual  probability  functions  [Stakgold  1998]. 

The  functional,  designated  by  (/,  0),  is  defined  by  the  integral  relationship 

(/,  f>)  =  J  f(x)  0  (x)  dx ,  (8) 

Rn 

which  may  be  multi-dimensional.  This  integral  spans  the  space  of  /  (x),  and  the  accessory 
functions  must  span  this  space  as  well.  The  derivative  of  the  functional  is  given  by 

=  (9) 


In  particular,  we  note  that  if  the  function  f  (x)  has  discrete  jumps  that  occur  at  points 
xp  then  its  functional  derivative  has  the  form 


d 

dx 


f 


—i 

dx 


+5>/,<s 


x 


Xi) 


(10) 


where  [f]  is  the  piecewise  continuous  derivative  and  A/,  is  the  signed  magnitude  of  the  ith 
jump.  We  shall  refer  to  this  as  the  jump  factor  or  function,  depending  on  form. 


Let  us  now  consider  the  differential  equation 


du 

dx 


(11) 


If  /  (x)  is  a  continuous  function,  then  we  may  find  a  classical  (or  strict)  solution  u  (x) 
that  obeys  the  differential  Equation  (11)  continuously  at  every  point  x.  If,  however,  /  is  a 
functional,  then  u  is  a  solution  of  the  differential  Equation  (11)  if  and  only  if 


(/>  <P) 


(12) 
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for  all  (f>  in  the  set  that  spans  the  space.  In  particular,  if  we  want  a  to  be  a  function, 
then  Equation  (12)  reduces  to  the  differential  Equation  (11)  in  the  weak  sense.  We  must 
note,  however,  that  any  weak  solution  that  has  m  continuous  derivatives,  where  rn  is  the 
dimension  of  the  space,  is  a  classical  solution.  For  our  purposes,  we  adopt  this  immediately 
as  an  assumption. 

What  this  translates  to  is  that  if  we  can  find  a  way  to  reduce  the  delta  functions  of 
Equation  (10)  to  having  a  weak  behavior,  then  we  reduce  that  equation  from  a  functional 
differential  equation  to  a  function  (or  classical)  differential  equation.  In  terms  of  our  search 
for  a  General  Theory  of  Rates,  this  differential  equation  is  a  (if  not  the)  fundamental  equa¬ 
tion. 


IV.  OVERVIEW  OF  LAPLACE  TRANSFORMS 


Our  second  body  of  theory,  which  is  particularly  useful  in  considering  deterministic  event 
times,  is  Laplace  Transforms  (LT)  [Doetsch  1961]. 

Assume  that  f  (t)  ,  a  function,  is  well  behaved  and  defined  for  t  >  0.  Then  we  may  define 
the  Laplace  Transform  of  /  (t)  as 


L  [/(f)]  =  / 

Jo 

=  /(«)• 

The  Inverse  Laplace  Transform  is  defined  as 


(13) 


-i  p+ifi 

f(t)  =  —  linr  /  estf  (s)  ds, 
2m  oo  J7_i/3 


(14) 


which,  unlike  the  Inverse  Fourier  Transform,  is  not  symmetric  and  is  considerably  more 
difficult  to  perform  than  the  transform  itself.  For  this  reason,  a  common  practice  with  LT 
is  to  invert  it  by  wimp  inspection  of  tables  of  transforms  rather  than  by  macho  integration. 


The  LT  also  obeys  the  convolution  theorem, 


fi  (t-t')  f2  ( t ')  dt 


Go 


Derivatives  have  the  transform 

dnf  ( t ) 
d,tn 


71—1 


=  S' 


“/  w-E 


=  fl  0)  f  2  (s)  • 


•  dn~i-lf(t) 


(15) 


j=0 


and  repeated  integrals  have  the  transform 


/•t  fti 

I  dtn..  /  dti/(G) 
Uo  Jo 


dtn~i~l 


=  s-nf  (s)  . 


(16) 


t= 0 


(17) 
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Also  useful  are  the  moments, 


dn  f  (s') 

L[tnf(t)]  =  (-ir^/.  as) 

The  expected  value  of  a  function  is  usually  given  as 

poo 

(t)=  t'f{t')dt\  (19) 

Jo 

assuming  f  (t)  is  defined  on  [0,  oo).  If  we  write  the  quantity  t'  f  (t/)  e~st'  dt\  which  is 
just  the  LT  of  t  f  (t),  then  we  see  immediately  that 

(t)  =  limL[tf(t)],  (20) 

and  by  Equation  (18),  we  may  make  the  association 

(()  =  -  lirn^hl.  (21) 

s^o  as 

While  multidimensional  LT  are  possible,  they  seldom  arise  so  we  shall  not  treat  with 
them  here. 

V.  PROBABILITY  THEORY 

If  an  event  is  stochastic,  it  has  a  Cumulative  Probability  Function  or  Cumulative 

Distribution  Function  (CDF)  that  we  designate  as  a  capital  Roman  letter  function  of 
time  (for  example,  P  (t))  and  possibly  other  parameters  or  variables.  This  introduces  an 
implicit  assumption  that  the  process  producing  the  event  is  continuous  and  mathematically 
well  behaved.  In  particular,  we  assume  that  the  CDF  is  at  least  m  times  differentiable. 
In  some  cases,  notably  with  the  combat  processes  in  the  next  chapter,  we  shall  deal  with 
discrete  probabilities,  but  they  are  still  assumed  to  be  well  behaved. 

Because  of  this  differentiability,  it  is  also  useful  to  deal  with  the  Probability  Density 
Function  (PDF)  that  we  designate  as  a  small  Roman  letter  function  of  time  (for  example, 
p(t))  and  possibly  other  parameters  or  variables.  The  CDF  is  the  probability  (we  use 
likelihood  as  a  synonym,)  that  an  observable  event  has  occurred  by  time  t,  while  the  PDF  is 
the  probability  per  unit  time  that  the  observable  event  has  occurred  within  the  infinitesimal 
interval  dt  of  time  t.  As  indicated  above,  these  are  related  as  derivative  and  integral,  so  that 
we  have  the  relationships, 

P(t)  =  (22) 

t 

P  (t)  =  J  p  (t')  dt' . 

o 


7 


In  addition  to  the  Guassian  and  NED  distributions  already  discussed,  we  also  have  the 
Gamma  Distribution  of  (integer)  order  n  (also  sometimes  called  the  nth  rank  Erlang  Distri¬ 
bution,)  given  by 

Pn  (t)  =  (atr1  e~at.  (23) 

[n  —  1)! 

This  distribution  has  expected  value  and  standard  deviation  of 

77 

(t)  =  (24) 

a 

w  —  > 
a 

which  is  sometimes  referred  to  as  a  Poisson  standard  deviation  since  this  distribution  occurs 
in  a  Poisson  process. 


Yet  another  distribution  we  want  to  consider  is  the  Lognormal  Distribution  given  by 


P(t) 


(In  (t)  -  ln(/i))2 
2  E2 


(25) 


which  is  just  a  normal  distribution  of  the  logarithm  of  the  random  variable  t.  This  distrib¬ 
ution  has  expected  value  and  standard  deviation 

(t)  =  fie-*2/2,  (26) 

a  =  n  e-s2  (e-s2  -  1), 

which  is  rather  complicated  and  makes  for  some  interesting  conversion  challenges.  Notably 
the  LT  of  the  lognormal  distribution  does  not  have  a  simple  representation.  It  is  easy  to 
calculate.  All  one  has  to  do  is  expand  the  exponential  e~st  in  an  infinite  series  about  zero  in 
either  variable.  The  integral  is  then  calculated  for  each  term  and  produces  an  exact  result. 
Unfortunately,  each  term  of  tn  produces  a  transform  term  of  e~n  *  /2.  Thus  the  resulting 
series,  which  is  the  LT,  is  well  defined  but  may  not  be  compactly  summed. 


For  deterministic  events,  we  shall  represent  the  PDF  of  the  event  by  the  Dirac  (strong) 
delta  function  defined  by  Equation  (1).  In  most  cases,  we  shall  be  dealing  with  the  translated 
delta  function  5  (t  —  t '),  which  has  non-zero  value  only  at  t  —  t! .  The  expected  value  and 
standard  deviation  of  this  translated  delta  function  distribution  are 


(t)  =  (27) 

a  —  0, 

which  are  what  we  should  expect  for  a  deterministic  event.  That  is,  the  event  has  a  definite 
time  of  occurrence. 


Two  elaborations  of  probability  theory  are  needed  to  describe  repeated  events  caused  by 
the  same  actor  and  events  caused  by  a  group  of  actors.  We  now  consider  each  in  turn. 


A.  Renewal  Theory 


Repetitive  actions,  that  is  repetitions  of  a  process,  result  in  a  series  of  repetitive  events. 
The  mathematical  formalism  used  to  describe  this  repetition  of  a  process  is  Renewal  The¬ 
ory  [Ross  1970].  ”A  renewal  (counting)  process  is  a  nonnegative  integer- valued  stochastic 
process  that  registers  the  successive  occurrence  of  an  event  during  a  time  interval,  where  the 
times  between  consecutive  events  are  positive,  independent,  identically  distributed  random 
variables”  [Taylor  Karlin  1998]. 

For  our  purposes,  we  are  primarily  concerned  with  the  repeated  nature  of  the  process  and 
shall  relax  the  identically  distributed  property  while  retaining  the  positive  and  independent 
properties.  This  generalization  is  consistent  with  both  the  mathematical  theory  of  renewals 
and  the  problem  at  hand.  In  the  absence,  as  yet,  of  detailed  knowledge  of  the  architecture 
of  the  rate  processes,  we  want  to  retain  as  much  generality  as  we  can. 

Renewal  Theory  counts  the  occurrence  of  repeated  events  whose  evolution  is  described 
in  terms  of  probability  functions,  usually  functions  of  time.  Each  event  is  referred  to  as  a 
renewal,  and  has  a  CDF  that  we  shall  designate  as  P^n>  (t),  where  the  parenthetical  super¬ 
script  indicates  the  number  of  the  renewal.  Events  are  sequential  -  that  is,  the  first  renewal 
occurs  before  the  second,  and  so  forth  -  and  independent.  The  total  number  of  renewals  may 
be  finite  or  infinite,  although  the  latter  situation  is  often  assumed  for  the  former  to  simplify 
the  mathematics.  We  adopt  this  convention  here,  thus  gaining  a  tacit  assumption  that  the 
ending  of  the  renewal  process,  for  any  of  the  reasons  above,  is  smooth  and  has  no  effect  on 
the  dynamics  of  the  process  representation. 

The  expected  number  of  renewals  by  time  t  is  expressed  by  the  Renewal  Function, 
which  is  defined  as 

OO 

P(R)  (t)  =  J^P{n)  (t),  (28) 

71—  1 

where  the  parenthetical  superscripted  “R”  indicates  renewal.  Because  the  CDF  are  assumed 
to  be  continuously  differentiable,  we  may  immediately  define  a  Renewal  Density  Function 
(RDF)  as 

OO 

/«>(*)=  £><”)(*).  (29) 

71=1 

Some  theorists  extend  the  summation  to  n  —  0  by  defining  a  zeroth  renewal  event  with  CDF 
P (t)  =  1,  from  which  we  may  deduce  a  PDF  for  this  renewal  of  p ®  (t)  —  5  (t),  where  5  (t) 
is  a  Dirac  delta  function. 

As  we  have  indicated,  renewal  theory  describes  repetitions  of  processes  culminating  in 
events,  counting  the  events  in  a  probabilistic  manner.1  These  repetitions  do  not  have  to  be 

Renewal  theorists  apparently  do  not  like  to  talk  about  events.  Since  our  interest  is  in  (at  least  stochas¬ 
tically)  observable  events  and  the  processes  that  produce  these  events  are  presumed  to  be  repeated  and  the 
mathematics  of  renewal  theory  describes  these  processes  and  events  well,  we  shall  take  the  pragmatic  but 
impure  approach  of  continuing  to  use  the  verbiage  of  renewal  theory. 
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identical  and  in  many  cases  are  not.  We  thus  define  the  PDF  of  the  nth  renewal  as  p(n)  (t). 
where  the  parenthetic  subscript  indicates  the  number  of  the  renewal.2  From  this,  we  may 
write  an  evolution  equation  of  a  renewal  in  terms  of  the  previous  one,  assuming  the  renewals 
are  Markovian  (that  is,  each  renewal  only  depends  on  the  previous  one,)  as 

t. 

P{n)  (; t )  =  f  G^{P^~l\  ..}  (t  -  t')  p(n)  (f)  dt'  (30) 

o 

where  G'^n-){..}is  the  initiator  of  the  initiation  of  the  renewal;  that  is,  it  represents  the  prob¬ 
ability  that  the  process  is  initiated  and  may  depend  on  the  number  of  the  renewal.  Notice 
that  this  integral  is  symmetric  in  interchange  of  the  time  dependence  of  the  two  components. 
This  will  play  an  important  role  in  subsequent  developments.  As  a  matter  of  convenience, 
we  shall  also  write  integrals  of  this  form  using  the  shorthand 

p(")  =  G<">*p(n).  (31) 

Two  special  situations  are  of  note.  First,  the  generator  may  be  simple  in  form,  depending 
only  on  the  preceding  renewal.  In  this  case,  Equation  (30)  has  the  form 

t 

P(n)  (t)  =  I  P(n"1}  (t  -  t')  p{n)  (t')  dt'.  (32) 

o 

We  shall  refer  to  this  as  a  simple  (or  linear)  renewal,  although  this  terminology  is  somewhat 
different  from  that  in  the  literature.  We  note  that  because  Equation  (32)  is  linear,  we  may 
immediately  write  the  evolution  equation  of  the  corresponding  renewal  PDF  as 

t 

P{n)  (t)  =  I p(n_1)  (t  -  if)  P(n)  (0  dt'.  (33) 

0 

Second,  if  all  of  the  renewal  PDFs  are  identical  {p(n)  (t)  =  p(t)),  then  Equation  (32) 
further  reduces  to 

t 

*■>(*)  =  /*■-*>  (t-Op(0<  (34) 

0 

with  corresponding  form  of  the  PDF.  We  shall  refer  to  this  as  a  simple,  uniform  renewal. 

Renewal  Theory  addresses  the  specific  concept  of  rate  through  Blackwell’s  Theorem.  This 
theorem  states  that 

lim  PW  (t  +  At)  -  (t)  —  — ,  (35) 

t— >oo  fl 

2 Even  in  General  Renewal  Theory,  all  of  the  (t)  are  identical  except  the  first.  What  we  are  doing 
here  is  a  generalization  of  the  mathematics  of  renewal  theory.  Our  calling  this  renewal  theory  is  loose,  but 
handy.  Our  mathematical  statement  of  renewal  theory  here  is  considerably  more  general  than  normally 
given.  In  most  cases,  our  application  of  this  theory  will  be  considerably  simpler  and  congruent  with  the 
usual  statement. 
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where:  /i  =  J0°°  p(t')  dt'  in  the  context  of  a  simple,  uniform  renewal.  The  quantity  1  / /i 

is  sometimes  called  the  rate  (of  the  renewal.)  Within  the  definition  of  the  renewal  density 
function,  this  reduces  to 

lim  ( t )  — *  — . 

t->oo  /i 

Recalling  that  Equation  (30)  is  symmetric  in  interchange  of  the  time  dependence  of  the 
two  functional  pieces  of  the  convolution,  we  may  also  write  a  differential  equation  form  of 
the  evolution  equations  as 

ftPMW  =  ..}(«)?(„,(  0)  (36) 

t 

-J  „}  (t  -  t’)  Tj ,(n|  ((-) 

o 

If  the  renewal  PDF  is  NED  with  rate  £,  then  Equation  (36)  further  reduces  to 

yP(n)  (t)  =  fG(n){P(T*_1), ..}  (t)  -  £P(n)  (t) .  (37) 

LLL 

This  is  an  exceedingly  useful  form. 

Two  further  considerations  are  worthy  of  note  here.  First,  the  combination  of  a  set  of 
renewals  is  itself  a  renewal  [Cox  Smith  1954] .  Second,  renewals  may  be  compounded. 


Remark  1.  Those  who  have  had  a  course  in  renewal  theory  will  likely  be  shaking  their  heads 
at  this  point.  This  is  not  the  standard  textbook  renewal  theory  for  a  couple  of  quite  good 
reasons.  First  of  all,  that  textbook  renewal  theory  isn’t  quite  general  enough  to  satisfy  our 
needs  for  some  of  the  problems  that  we  want  to  solve.  Second,  the  constraints  that  we  need 
to  apply  are  somewhat  different  from  those  normally  applied  in  the  textbook  presentations. 
Yes,  this  is  very  different  from  vanilla  renewal  theory.  It’s  different  enough  that  I  might  have 
called  it  repetition  theory  except  that  then  someone  would  have  accused  me  of  plagiarism. 


B.  Order  Statistics 

By  the  postulates  of  General  Rate  Theory,  we  may  deal  with  events  that  are  caused  by  a 
collection  of  agents.  The  fact  that  we  have  a  collection  may  mean  that  these  agents  interact 
in  some  fashion.  In  particular,  we  know  that  this  collectivity  has  an  effect  on  the  probability 
of  the  event  even  when  the  agents  act  independently.  The  means  for  considering  this  is  Order 
Statistics. 

Where  we  used  Basic  Probability  Theory  and  its  extension  Renewal  Theory  to  describe 
the  dynamics  of  events  caused  by  a  single  element,  we  use  Order  Statistics  to  describe  the 
behavior  of  the  several  as  a  group.  The  basic  idea  of  Order  Statistics  is  quite  simple,  but  its 
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implementation  can  be  somewhat  difficult,  especially  in  terms  of  the  basic  calculations  and 
computations  that  need  to  be  performed.  For  this  reason,  it  is  useful  to  consider  pictures. 

The  notion  of  pictures,  as  we  use  it  here,  is  taken  from  physics,  specifically  quantum 
mechanics.  In  quantum  mechanics,  there  are  two  pictures  or  interpretations  of  how  the  time 
dependence  of  the  mathematical  representation  is  viewed.  The  composition  of  quantum 
mechanics  and  indeed  the  nature  of  these  pictures  is  not  directly  relevant  here  except  for  the 
fact  that  the  two  pictures  are,  within  the  restrictions  of  their  formulations,  equivalent  and 
give  identical  observable  results. 

We  want  to  consider  two  pictures  of  probability,  which  we  may  call  the  implicit  and 
explicit  temporal  representation  pictures.  The  explicit  picture  is  the  more  familiar.  In  this 
picture,  the  probability  functions  are  seen  as  dependent  variables  of  (nominally)  time.  In  the 
explicit  time  picture,  we  operate  directly  on  the  time  variable.  In  the  implicit  picture,  we 
operate  directly  on  the  probability  functions  as  variables.  This  requires  that  the  probability 
functions  (at  least  the  CDFs)  are  monotonically  increasing  over  time  and  that  there  is  a 
one-to-one  correspondence  between  probability  and  time.  This  is  not  particularly  difficult, 
but  there  is  one  idea  that  we  must  accept.  In  this  picture,  we  may  speak  of  the  idea  that  a 
probability  may  have  moments  (for  example,  expected  value)  and  because  of  the  one-to-one 
relationship,  this  implies  a  time. 

In  particular,  this  also  means  that  we  must  assume  that  the  probability  functions  are 
invertible.  That  is,  if  P  (t)  is  the  CDF  that  some  (observable)  event  has  occurred,  then  there 
exists  an  inversion  t  (P)  that  is  the  unique  time  (not  necessarily  finite)  when  that  probability 
value  occurs. 

We  may  now  take  up  our  consideration  of  Order  Statistics. 

Consider  a  collection  of  N  elements,  each  engaged  in  the  same  stochastic  process  (simple 
or  complex)  that  results  in  an  observable  event.  We  thus  expect  to  observe  N  events,  if  we 
observe  for  a  sufficiently  long  time.  Assuming  we  know  the  stochastic  mechanics  of  individual 
elements,  given  by  a  CDF  P  (t).  then  we  naturally  ask  the  question  of  when  do  the  events 
occur.  (Or  equivalently,  What  does  the  timing  of  a  few  events  tell  us  about  the  collection?) 
The  way  we  address  this  is  by  adopting  the  implicit  time  picture  and  back  into  the  problem. 

Let  us  specifically  ask  the  question,  What  is  the  probability  that  the  Prevent  has  oc¬ 
curred?  If  we  write  the  unit  function 

Z  =  [P+(1-P)]N  (38) 

=  1, 

then  a  simple  application  of  the  binomial  theorem  allows  us  to  write  this  function  as  a  series 

N  /  /y\ 

Z  =  (  )  (!  -  p)  pr ,  (39) 

r=0  '  ' 
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whose  individual  terms  are  the  probabilities,  within  the  collection  environment,  that  the 
rtheve nt  has  occurred.  Thus,  for  explicitness, 

Pr:N=  ^  (1  ~P)N~rPr,  (40) 

is  the  probability  that  the  rtheve nt  out  of  N  has  occurred.  Because  the  component  proba¬ 
bilities  are  time  dependent,  this  probability  is  also  time  dependent. 

Since  the  CDF  exists  and  is  time  dependent,  this  implies  that  its  PDF,  its  first  derivative 
with  respect  to  time,  exists  as  well.  Taking  this  derivative  is  not  completely  straightforward. 
Because  the  term  (1  —  P)N  '  represents  the  probability  that  the  event  has  not  occurred  for 
N  —  r  of  the  elements,  it  does  not  contribute  to  the  derivative.  Thus,  the  PDF  of  the  rth event 
is  just 

Pr:N=  (^)r  (1  ~  Pf~r  P'-'p,  (41) 

where  p  =  4P.  This  is  admittedly  a  very  terse  presentation  and  the  student  is  referred  to 
[David  1970]  for  details. 

As  before,  this  PDF  is  time  dependent,  and  we  may  consider  either  the  expected  time 
the  event  occurs  in  the  usual  explicit  time  picture  or  the  expected  value  of  the  individual 
element  CDF  in  the  implicit  time  picture.  This  demonstrates  the  value  of  the  two  pictures. 
The  expected  time  of  the  event  is 

OO 

tr:N  =  J  t  Pr-.N  (t)  dt  (42) 

0 

OO 

=  f  t(l-P  ( t))N~r  rP  (t)r_1  p  ( t )  dt, 

o 

which,  given  the  explicit  nature  of  the  probability  functions,  is  likely  to  be  quite  complicated. 
The  expected  probability  of  the  event  is 

i 

Pr:N  =  J  P  dPr:N  (P)  (43) 

0 

1 

=  /  P  (i P)N~r  rPr~1dP 

o 

r 

N  +  V 

which  is  considerably  simpler.  Of  course,  we  must  still  compute  the  value  of  time,  but  this 
is  given  by  the  assumption  of  inversion  of  the  CDF,  t  (yryi)  •  Note  that  we  have  made  use 

of  the  differential  relationship  pdt  =  dP  and  the  integral  expression  x (1  —  x) vl  dx  = 

rtdrp) 

T(n+i>) 
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VI.  ON  TO  THE  BASICS 


The  previous  sections  have  been  pretty  heavy  stuff,  so  let’s  cut  through  some  of  the 
fog.  What  we  are  all  about  is  finding  means  for  representing  the  counting  of  changes  in 
populations  of  discrete  events  in  a  continuous  manner.  Our  starting  point  for  this  is  Equation 
(10), 

L  -1  l 

which  we  shall  refer  to  as  the  basic  jump  differential  equation.  In  most  cases,  we  are  going 
to  be  dealing  with  situations  where  we  start  with  discrete  events  and  want  to  end  up  with  a 
continuous  representation  of  these  events  (as  shown  in  Figure  1.)  As  a  result,  in  general  the 
continuous  derivative  term  of  this  equation  will  be  zero.  Therefore,  we  are  going  to  discard 
this  term  and  reduce  this  equation  to 

=  '52&fiS(x-xi).  (44) 

i 

We  note  that  even  should  we  have  continuous  derivative  contributions,  it  is  simple  to  put 
them  back  in,  so  we  give  up  this  modicum  of  generality  with  the  understanding  that  we  may 
re-insert  it  at  any  time,  and  by  seizing  this  speciality,  we  gather  some  degree  of  simplicity 
and  thereby  clarity. 

As  a  matter  of  course,  we  shall  refer  to  the  specific  events  as  jumps  because  we  presume 
them  to  be  impulsive.  These  jumps  are  represented  by  the  terms  of  the  series  on  the  right  side 
of  the  above  equations.  We  shall  refer  to  the  magnitudes  of  these  jumps,  A/*,  as  jump  factors 
or  jump  functions.  The  positions  of  the  jumps  or  events  are  the  Xi,  which  we  commonly  take 
to  be  times,  but  which  may  be  spatial  instead.  Regardless,  these  jump  positions  or  times 
are  taken  to  be  fundamentally  independent  of  the  jump  differential  equation  although  they 
may  result  from  some  phenomena  that  are  described  by  some  conjugate  body  of  theory  or 
experiment.  To  further  clarify  this,  we  rewrite  this  discrete  event  jump  differential  equation 
as 

=  (45) 

% 

to  explicitly  show  that  we  shall  generally  consider  the  jump  positions  to  be  times. 

It  is  somewhat  useful  in  developing  this  theory  of  rates  to  introduce  a  bit  of  taxonomy 
based  on  the  nature  of  these  jump  factors/functions  and  times/positions.  This  taxonomy  is 
summarized  in  the  Table  1. 

The  taxonomy  differentiates  the  situations  of  Deterministic  and  Stochastic  behavior  of 
the  jump  functions  and  times  by  numbered  type.  Thus,  for  example,  a  type  one  jump 
differential  equation  is  characterized  by  both  the  jump  function  and  the  jump  times  being 
deterministic,  while  a  jump  differential  equation  of  type  four  is  characterized  by  both  jump 
functions  and  jump  times  being  stochastic.  As  we  shall  see,  this  taxonomy  is  presented 
largely  to  provide  a  framework  for  presenting  methods  for  building  the  continuous  forms  of 
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the  jump  differential  equations  from  the  discrete  forms.  The  distinction  being  arbitrary,  we 
are  erecting  no  barriers  among  differential  equations  that  may  be  of  mixed  type. 


Table  1.  Jump  Taxonomy 


Type 

Jump  Function  Form 

Jump  Time  Form 

1 

Deterministic 

Deterministic 

2 

Stochastic 

Deterministic 

3 

Deterministic 

Stochastic 

4 

Stochastic 

Stochastic 

Our  goal,  as  we  have  indicated,  is  to  describe  how  we  may  start  with  discrete  event 
differential  equations,  such  as  the  above,  and  transform  them  into  a  continuous  form.  This 
taxonomy  gives  us  a  means  of  presenting  different  methods  of  transformation. 

As  a  rule,  in  this  presentation,  we  shall  assume  that  the  jump  functions  do  not  depend 
on  the  index  of  the  jump;  that  is,  they  either  depend  only  on  time  or  are  independent  of  the 
jump  index.  As  a  result,  we  may  usually  separate  the  jump  function  from  the  series. 

VII.  TRANSFORMING  THE  DELTA  FUNCTIONS 

Let  us  now  return  to  Equation  (45),  and  as  we  have  indicated,  presume  that  the  jump 
functions  do  not  depend  on  the  index  of  summation.  In  this  case,  we  may  rewrite  the  right 
side  of  the  equation  as 


Y  Afi  S(t-  U )  =  A f  (t)  Y^S(t-  U )  ,  (46) 

i  i 

where  we  have  assumed  that  the  jump  function  is  a  function  of  time.  As  we  shall  see  shortly, 
this  does  not  diminish  the  overall  utility  of  the  theory. 

As  we  have  also  indicated,  we  want  to  find  a  means  of  transforming  this  discrete  summa¬ 
tion  over  events  into  some  continuous  representation.  We  now  proceed  to  examine  exactly 
this.  Our  discussion,  for  simplicity  and  clarity,  is  limited  to  a  single  or  sequence  of  event 
times.  Combinations  of  sets  or  sequences  complicate  the  analysis  but  are  ready  extensions 
of  the  theory. 

A.  Equally  Spaced  Deterministic  Event  Times 

We  start  our  presentation  by  considering  a  set  of  event  times  that  are  infinite  in  number 
and  equally  spaced.  As  a  result,  we  may  rewrite  Equation  (46)  explicitly  as 

OO 

a/  (t)  *0 =  (*)  >  (47) 

i  i= 0 
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where  we  have  arbitrarily  assumed  that  the  jumps  begin  at  index  zero.  If  we  define  =  i  At, 
to  assure  that  the  events  are  evenly  spaced  in  time,  then  this  becomes 

OO 

Af  (*)  ti)  =  Af  (t)  ~ i  At)  •  (4§) 

i  i=0 


Let  us  now  take  the  Laplace  Transform  of  this  equation 


L 


Af  {t)^2S(t-i  At) 

i= 0 


i= 0 


e~stAf  (t) 


5  (t  —  i  A t)dt, 
8  (t  —  i  A t)dt, 


XV“  A‘Af  (’ A() . 

i= 0 


(49) 


where  we  have  had  to  pull  the  jump  function  back  into  the  summation  because  it  now 
depends  on  event  number  courtesy  of  the  form  of  the  delta  functions.  Now,  if  the  jump 
function  A /  (t)  is  strictly  polynomial  (and  preferably  of  low,  positive  order,)  then  this  series 
is  summable,  although  for  even  reasonably  small  powers  of  i,  this  summation  will  be  difficult. 
Accordingly,  we  would  prefer  not  to  take  this  approach  to  making  the  transformation  from 
discrete  to  continuous  representation. 


Accordingly,  we  limit  ourselves  to  taking  the  Laplace  Transform  of  just  the  series  of  delta 
functions 


8  (t  —  i  At) 


i=0 


=  E« 

2—0 


— si  At 


(50) 


which  we  immediately  see  is  summable  as 


L 


y  5  (t  —  i  At) 

_t=o 


1 

1  -  e~sAt  ’ 


(51) 


since  Equation  (50)  is  a  simple  geometric  series  [Gradshteyn  Ryzhik  1980,  Section  1.112 
Equation  1,  p.  21.]. 


As  it  stands,  this  result  is  not  an  improvement  and  has  gotten  us  no  closer  to  a  continuous 
representation.  If  however,  we  expand  the  exponential  to  first  order, 


e~sAt  ~  1  -  sAt  +  HOT, 


then  we  may  approximate  Equation  (51)  as 


L 


y  5  (t  —  i  At) 

_i=0 


sAt’ 


(52) 


(53) 
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which  we  may  invert  by  inspection  and  obtain 

OO 

8  (t  ~  i  At)  ~  — ,  (54) 

i=0 

which  is  equivalent  to  an  infinite  time  limit  average  (ala  Blackwell’s  Theorem  as  we  shall  use 
subsequently  in  considering  Type  3/4  jump  differential  equations.)  That  is,  we  spread  the 
magnitude  of  the  jumps  over  the  entire  duration  of  jumps. 

This  result  transforms  the  jump  differential  equation  to  a  continuous  form, 

(55) 

d,tJ  At  v  ’ 

While  this  is  a  very  useful  result,  it  is  limited  in  that  we  cannot  extend  the  method  to 
higher  order  approximations  of  the  exponential  without  risking  instability.  Accordingly,  let 
us  now  consider  that  instead  of  the  first  jump  occurring  at  the  origin,  that  it  occurs  after  an 
increment  of  time.  Thus,  we  rewrite  Equation  (47)  as 

OO 

a/  (t)  s  (t  -  u )  =  a  /  (t)  (56) 


We  now  follow  the  same  method  as  before,  assuming  equal  increment  times,  and  obtaining 


L  8  (t  —  i  At) 


1  -  e~sAt 

e-sAt 

~  1  —  e~sAt  ’ 

1 

“  esAt  _  i  ■ 

This  may  be  expanded  as  before.  A  first-order  expansion  gives  us  the  same  result  as  before, 
Equation  (53).  Unlike  that  expansion,  which  was  limited  to  first  order,  we  may  expand 
Equation  (57)  to  higher  order.  Thus,  an  expansion  to  second  order  gives  us 


L  8  (t  —  i  At) 


sAt-\- 


s2  At2  ' 


This  can  be  separated  using  a  partial  fraction  expansion,  and  inverted  by  inspection,  giving 


Y^8{t-iAt)  [1  -  e  2t/At\  . 


We  immediately  note  that  in  the  limit  t  — »  oo,  this  result  becomes  identical  to  Equation  54. 
Similarly,  it  gives  us  a  continuous  representation  rate  differential  equation 
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As  a  means  of  comparison,  we  plot  the  rate  functions  arising  from  the  first  three  order 
expansions  of  the  LT  of  the  delta  function  series,  Equation  (57),  in  Figure  2.  The  third- 
order  expansion  is  actually  negative  out  to  a  time  of  about  0.7,  but  this  is  clipped  in  the 
plot.  It  may  readily  be  seen  that  the  detail  of  representation  of  the  delta  functions  increases 
with  expansion  order.  In  a  similar  manner,  we  present  the  integrals  of  the  first  two  order 
expansions  versus  the  step  trajectory  in  Figure  3.  We  may  readily  see  that  the  type  of 
“smoothing"  we  displayed  earlier  in  Figure  1  occurs  here,  but  we  are  able  to  see  greater 
detail  of  the  effects  of  expansion  of  the  LT  of  the  event  delta  function  series.  The  first-order 
expansion  provides  a  smooth  representation  that  hits  the  vertical  risers  of  the  “steps,"  while 
the  second  order  expansion  hits  at  the  mid  points  of  the  risers. 
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Figure  2.  Comparisons  of  First  Three  Order  Expansions  of  Series  of  Event  Delta  Functions 


Before  proceeding,  it  seems  worthwhile  to  make  a  comparison  of  this  method  of  devel¬ 
oping  the  continuous  representation  by  expansion  of  the  Laplace  Transform  of  the  discrete 
representation  with  a  more  direct  method.  The  direct  method  we  shall  consider  is  based  on 
the  expression  of  delta  functions  as  Gaussians,  Equation  (4).  As  already  noted,  the  Gaussian 
is  a  weak  delta  function,  which  means  that  in  the  limit  a  — *  0,  the  Gaussian  acts  like  a  strong 
(or  Dirac)  delta  function.  Our  approach  here  is  the  opposite  of  this.  Instead  of  making  a  get 
smaller,  we  let  it  get  larger.  This  is  shown  in  Figure  4.  As  sigma  is  allowed  to  increase,  the 
discrete  Gaussians  coalesce.  Thus,  we  obtain  the  same  type  of  continuous  representation. 
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Figure  3.  Comparisons  of  Integrals  of  First  Two  Order  Expansions  of  Series  of  Event  Delta 
Functions  With  Actual  Step  Trajectory 


B.  Unequally  Spaced  Deterministic  Event  Times 

The  problem  we  encounter  in  dealing  with  unequally  spaced3  event  times  is  that  the 
transformation  techniques  used  above  to  go  from  discrete  to  continuous  representation  are 
not  generally  applicable.  To  illustrate  this,  we  return  to  Equation  (56) 

OO 

Af  (f)  &  (t  ~  U)  =  A /  (t)  ^ps(t-ii)  , 

i  2=1 

where  we  have  assumed  the  jump  function  to  be  independent  of  index. 


is 


As  a  result,  we  need  only  consider  the  series  of  delta  functions,  whose  Laplace  Transform 


L 


_  2=1 


X>“*- 

2=1 


(61) 


By  assumption  there  is  no  simple  linear  relationship  between  index  and  event  times,  and 
we  cannot  simply  sum  this  series  as  we  did  before.  If  we  examine  the  event  times  and 


'Deterministic.  By  their  nature,  stochastic  event  times  are  unevenly  spaced  but  are  treated  in  a  different 
manner. 
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Figure  4 ■  Effect  of  Increasing  Sigma 


find  that  the  difference  in  adjoining  event  times  is  almost  constant,  then  we  may  be  able 
to  approximate  them  as  evenly  spaced,  or  if  the  difference  is  almost  constant  except  for  a 
disordered  but  small  amount,  we  might  be  able  to  approximate  them  as  evenly  spaced  with 
a  small  stochastic  component  that  can  the  transformed  using  the  techniques  described  in  the 
next  subsection. 

Generally  however,  we  cannot  sum  Equation  (61)  in  a  simple  manner  so  we  turn  rather 
quickly  to  transforming  the  series  term  by  term  using  the  techniques  described  in  the  previous 
subsection.  Thus, 


E^-E 


E— E 

1  +  sti 

1=1 

oo 


(62) 
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which  we  may  invert  as 


OO 


OO 


L~' 


y  -si. 


(63) 


Although  elegant,  this  suffers  from  being  tied  to  the  origin.  As  a  result,  density  of  rate 
is  often  too  bunched  at  short  times.  This  is  illustrated  in  Figure  5  where  we  compare  this 
approximation  with  a  set  of  uneven  event  times.  The  event  times  in  this  case  are  generated 
from  the  deciles  of  a  Gaussian  with  mean  and  standard  deviation  selected  to  assure  all  the 
event  times  are  nonnegative.  The  different  heights  of  the  delta  functions  are  an  artifact  of 
the  small  standard  deviation  (a  =  0.01)  and  the  calculation  increment.  The  general  failure 
of  the  approximation  to  capture  the  behavior  of  the  events  is  obvious  in  this  case. 


Figure  5.  Comparison  of  Uneven  Event  Times  Delta  Functions  And  Term-By-Term  LT 
Expansion 


As  a  result,  we  are  often  forced  on  the  more  difficult  task  of  taking  some  representation 
of  the  delta  function,  relaxing  its  strength  (increasing  its  weakness),  and  looking  for  some 
pattern  of  behavior  that  can  be  exploited  in  forming  a  continuous  (approximate)  represen¬ 
tation.  For  example,  if  we  take  the  same  event  times  as  used  in  Figure  5,  represent  the  delta 
functions  as  Gaussians  (which  we  did  in  that  figure  for  plotting  purposes  anyway,)  and  allow 
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the  standard  deviation  (er)  to  increase,  then  we  find  a  pattern  coalescing.  This  is  shown  in 
Figure  6. 
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Figure  6.  Comparison  of  Uneven  Event  Times  as  Gaussians  of  Different  Standard  Deviations 


We  first  note  that  to  avoid  hiding  details  because  of  maximum  amplitude  changes,  we 
have  plotted  the  same  representation  of  the  delta  functions  of  Figure  5  on  one  vertical  axis 
and  the  other  two  representations  on  a  second  vertical  axis.  The  amplitude  variation  of  the 
small  standard  deviation  representation  is  again  an  artifact  of  calculation.  When  we  relax 
the  standard  deviation  from  0.01  to  0.1,  then  we  see  the  emergence  of  an  envelope  with  the 
peaks  sticking  through.  If  we  further  relax  the  standard  deviation  to  0.3,  then  the  curve 
smooths  and  we  recognize  it  as  a  Gaussian,  which  is  what  we  should  expect  given  the  event 
times  being  the  deciles  of  a  Gaussian. 

C.  Stochastic  Event  Times 

As  indicated,  the  times  of  deterministic  events  are  represented  by  strong  (Dirac)  delta 
functions  and  we  form  the  continuous  representation  by  weakening  them  in  some  manner. 
These  need  no  further  consideration  at  this  point.  The  behavior  of  the  stochastic  events 
does  require  development,  however.  Being  stochastic,  we  intuit  that  the  event  times  are 
represented  by  delta  functions  that  are  inherently  weak  in  some  sense.  To  examine  this,  we 
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rewrite  our  jump  differential  equation,  Equation  (44),  with  probability  as  the  independent 
variable.  Thus, 

df  (t)  =  fi5(P-  Pi)  dP  0 1 ) .  (64) 

i 

This  telegraphs  that  we  are  going  to  use  the  implicit  temporal  representation  picture,  at 
least  initially. 

Now  consider  a  collection  of  N  elements.  These  elements  cause  events  that  define  the 
solution  /  (t)  of  the  differential  equation.  We  start  with  the  restriction  that  each  element 
causes  (at  most)  one  event.  Under  this  restriction,  basic  probability  theory  defines  the  CDF- 
PDF  of  the  event  caused  by  an  individual  element.  Thus,  P  (t)  represents  the  cumulative 
probability  that  the  event  caused  by  an  individual  element  has  occurred. 

Because  the  elements  are  parts  of  a  collective,  the  individual  events  are  described  by 
order  statistics.  As  a  result,  the  PDF  of  the  rth  of  N  events  is  given  by  Equation  (41),  which 
we  now  write,  with  differential,  as 


pr-.N  (t)dt  —  (  J  (1  —  P)N  r  rPr  1dP.  (65) 

It  may  be  noted  that  the  left  side  is  explicitly  a  function  of  time  and  the  right  side  is 
explicitly  a  function  of  individual  element  probability.  The  right  side  may  be  considered 
to  be  the  PDF  as  a  function  of  probability  (times  a  differential).  We  concentrate  on  the 
non-differential  portion.  This  is 


Pr-.N  (. P )  =  (N)  (1  -  P)N  r  rPr~1 .  (66) 

Let  us  now  examine  this  function  of  the  individual  element  CDF.  We  plot  selected  functions 
in  Figure  7.  These  functions  clearly  have  the  appearance  of  weak  delta  functions  in  the  sense 
of  our  earlier  discussion,  but  this  needs  to  be  demonstrated  in  the  context  of  Equation  (5). 
To  do  this,  we  return  to  Equation  (66),  and  rewrite  it  somewhat  as 

Pr-.N  ( t )  d.t  =  (1  -  p)N~r+1  rPTd  In  ^  ^  p)  '  (67) 


We  select  the  variable  y  =  In  (p^-p)  f°r  two  reasons.  First,  it  gives  us  natural  limits  of  ±oo, 
and  second,  as  we  shall  see,  it  reproduces  the  expected  values  of  Equation  (41).  Removing 
the  differential,  we  write  the  quantity  as 


Pr-.N 


1  -  P 


re(JV-r+l)  ln(l-P)+r  ln(P) 


(68) 
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Since  we  know  these  functions,  at  least  in  the  P  picture,  are  peaked,  we  may  expand  the 
exponent  in  a  McWheeny  or  (pseudo)  Saddle  Point  Expansion.  This  is  just  a  Taylor  series 
expansion  where  the  expansion  point  is  the  extremum,  normally  the  maximum.  By  this,  we 
mean  that  we  expand  the  exponent  in  y  =  In  (jz-p)  and  set  the  first  derivative  to  zero.  This 
determines  the  expansion  point.  For  simplicity  of  presentation,  we  are  only  going  to  keep 
the  first  two  derivative  terms,  recognizing  that  this  results  in  some  error,  but  we  only  want 
to  demonstrate  the  weak  delta  function  behavior  of  this  class  of  functions. 


Figure  7.  Order  Statistics  PDF  Examples 


We  now  write  out  the  expansion  explicitly  to  second  order, 


F 

d 

dyF 


<P_ 

dy2^ 


(N  —  r  +  1)  In  (1  —  P)  +  r  In  (P) , 
d 

<9  hi  (prp)  ^ 
dP  d 
din  (y^p)  dPr' 


P(l-P) 


(N  —  r  +  1) 


—  1  r 
1  -  P  +  P 


—  (N  +  1)  P  +  r; 

d  f  d  \  _  dp  d 
\dyF )  din  (j^p)  dP 


-  (N  +  1)  P  (1  -  P) . 


(69) 
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The  third  and  fourth  derivatives  are 

=  —  (N  +  1)  P  (1  —  P)  (1  —  2P) ,  (70) 

~^F  =  -(IV  +  l)P(l-P)(l-6P  +  6P2). 

The  pattern  is  not  immediately  obvious,  but  at  some  point  the  derivative  changes  sign 
(depending  on  the  value  of  P)  and  thereafter  becomes  oscillatory.  This  series  thus  also 
becomes  uniformly  convergent  and  assures  that  the  behavior  we  are  capturing  in  this  two- 
term  expansion  is  trustworthy. 


If  we  now  require  the  first  derivative  to  be  zero,  we  obtain  the  expansion  point  P  = 
which  is  just  the  expected  value  of  the  probability,  Equation  (43).  Substitute  this  into 
Equation  (70), 


F 


d_ 

dyF 

<P_ 
dy2  F 


(N  —  r  +  1)  In  (  1  — 
(N  —  r  +  1)  In 

T 

~(N  +  1) 


N  +  1 
N  —  r  +  1 
N  +  l 


+  r  In 
+  r  In 


r 


N  +  l 
r 

N  +  l 


0; 

-  (JV  +  1) 


JV  +  1 

r  N  —  r  +  1 


N  +  l  N  +  l  ’ 
r  (N  —  r  +  1) 


N  +  l 


which  gives  us  the  expansion 


F(y) 


~  (N  —  r  +  1)  In 


N  —  r  +  1 


N  +  l 

r  (N  —  r  +  1)  [y-Wiv^Tl)] 
N  +  l  2 


+  r  In 
2 


N  +  l 
+  O  (!/S)  . 


(71) 


(72) 


If  we  now  define  y  =  In  ( JV_rr+1)  ,  o2 


r(N-r+ 1)  ’  then  we  may  rewrite  Equation  (72)  as 


/  (x)  ~  (N  —  r  +  1)  In 


/  IV  —  r  +  1  \ 

V  N  +  l  ) 


+  r  In 


N  +  l 


(: y-yf 

2a2 


(73) 


and  thus  rewrite  Equation  (68)  as 


Pv.n  (y) 


^/Y,(  V  r  •  i )  in(  vv !  i 1 )  ■  r  i"(  v,|  i ) 


(74) 
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This  equation  may  be  more  explicitly  rewritten  and  rearranged  as 


Pr:N  ( y )  = 


N\  r  (N  —  r  +  1) 


r!  (N  —  r)!  (N  +  1)N+1 
To  make  factors  cancel  correctly,  we  further  rewrite  this  as 


N— T+l  r  2 

T  ( y-y ) 

P  2cr2 


(IV  +  1)!  r  (N  -  r  +  1)  (TV  —  r  +  l)iV~'r+1  +  (y-g)2 


Pr:N  (y)  ^ 


(IV  +  1)  r!  (AT  -  r  +  1)!  (N  +  1) 
and  make  use  of  Stirling’s  Approximation  for  the  factorial 


N+l 


2cr2 


(75) 


(76) 


n 


!  ~  nn+1/2e-"V^. 


(77) 


Now  consider  just  the  terms  on  the  right  side  of  Equation  (76)  leading  the  exponential, 

(N  +  l)!  r  ( N-r  +  1)  (IV  —  r  +  l)^-^1  rr 
(N  +  1)  r!  (N  -  r  +  1)!  (iV  +  l)^"1"1 

(IV  +  i)Ar+1+1/2  e~N~1r  (N  -  r  +  1)  (IV  -  r  +  l)^-^1  rr 

V^7r  (IV  +  1)  rr+l/2e~r  (N  —  r  +  i)N~r+1+1/‘1  e-N+r-i  (N  +  i)^1 
(IV  +  1)1/2  e~N~1r  (AT  —  r  +  1) 

V^7r  (IV  +  1)  r1/2e_r  (.V  +  r-  1)1/2  e-^”1 
(IV  +  1)1/2  r  (IV  —  r  +  1) 

V^7r  (N  +  1)  r1/2  (IV  —  r  +  1)1/2 
1  lr  (N  —  r  +  1) 

V/27r  V  IV  +  1 
1 


This  reduces  Equation  (76)  to 


(2/)  ^ 


1  (y-y)2 


(79) 


which  we  immediately  recognize  as  a  Gaussian.  The  standard  deviation  a  0  (  Jn  )  so 
that  cr  — >  0  weakly  as  N  — >  oo.  This  satishes  the  definition  of  a  weak  delta  function,  and  we 
may  thus  rewrite  Equation  (68)  as 


Pr-.N  (t)  (lt  =  6(y  —  y)  dy, 


(80) 


in  this  weak  sense.  The  relationship  of  the  Saddle  Point  Expansion  to  the  original  functions 
is  shown  in  Figure  8.  The  agreement  is  inexact  because  we  truncated  the  expansion.  It  is 
also  worthwhile  to  note  that  this  expansion  is  valid  for  N  >  0.  This  is  important  since  it 
allows  us  to  address  collections  of  cause  agents  of  population  one. 
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Figure  8.  Order  Statistics  PDFs  and  Their  Saddle  Point  Expansions 


VIII.  RATE  DIFFERENTIAL  EQUATIONS 

We  now  complete  the  derivation  of  the  stochastic  rate  differential  equations,  review  the 
deterministic  rate  differential  equations  for  evenly  space  event  times,  and  compare  the  two. 

A.  Stochastic  Rate  Differential  Equations 

We  now  combine  this  result  with  our  jump  differential  Equation  (64), 


(81) 


N 


The  differential  equation  is  now  a  function  of  time,  and  the  summation  over  jumps  has  be¬ 
come  explicit.  This  is  the  basic  stochastic  rate  differential  equation  (except  for  the  piecewise 
continuous  derivative  term.)  Although,  we  do  not,  at  this  point  in  the  development,  know 
the  explicit  form  of  the  jump  function,  some  further  elaborations  are  still  possible. 

When  the  jump  function  is  not  explicitly  a  function  of  i.  we  may  pull  the  jump  function 
out  of  the  summation.  Further,  if  we  substitute  the  exact  expression  for  the  order  statistics 
PDF,  Equation  (65),  we  obtain 


r—1 


(82) 
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(83) 


Writing  out  the  binomial  coefficient  and  canceling  the  r  term  yields 


TV 

df  (; t )  =  A  f  (t)  Tf 

r=l 


A! 

W 


—  (l-P)JV_rPr"1pdt 
r)! 


We  now  make  the  change  of  summation  index  s  —  r  —  1,  so  that  Equation  (83)  becomes 

TV— 1 


dfit)  =  A 


N\ 


1)! 

TV-1 


(1  ~P) 


TV— s— 1 


Pspdt 


(84) 


^S!(JV-S-1)! 


a/  w  .,AV  TT,  (i  -  pf—'  P ’  dt 


=  A  f(t)Npdt. 


This  collapse  results  from  recognition  that  the  middle  representation  of  Equation  (84)  is 
itself  a  Binomial  Expansion.  If  we  now  recognize,  based  on  our  earlier  development  of  the 
real-valued  binomial  expansion,  that  N  may  be  time  dependent  and  that  the  individual 
element  PDF  is  time  dependent,  we  may  write  this  as 

jf  f  (t)  =  A/  (t)  p  (t)  N  (t) .  (85) 

This  is  a  general  result  for  the  case  when  the  jump  function  does  not  depend  on  the  event 
number. 


Before  concluding  this  section,  two  other  cases  need  to  be  considered.  The  first  is  when 
the  jump  process  is  repeated  and  described  by  renewal  theory.  In  this  case,  we  make  use  of 
the  independence  of  the  individual  renewals  and  extend  it  to  the  order  statistics.  Equation 
(81)  becomes 

OO  N 

df  (t)  =  Afrn)  (t)  dt,  (86) 

n= 1  r— 1 

where  A  and  p^N  (t)  are  now  the  jump  function  and  order  statistic  PDF,  respectively  of 
the  rth  event  of  N  of  the  nth  renewal  or  repetition.  If  the  jump  function  does  not  depend  on 
event  number,  but  does  depend  on  renewal  number,  then  this  equation  may  be  written  as 

OO  N 

df  (t)  =  Af(n)  Pr?N  (*)  dt-  (87) 

n= 1  r= 1 

This  equation  collapses  as  before  to  give 

OO 

df  (; t )  =  A  f{n)p(n)N  dt.  (88) 

n=  1 

If  the  jump  function  depends  on  neither  event  nor  renewal  number,  then  this  becomes 

OO 

d/(i)  =  A/(f)5><  n^Ndt,  (89) 

n=  1 
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which  reduces  to 


df[t)  =  Af{t)p{R)Ndt.  (90) 

This  is  the  form  that  we  shall  use  subsequently  for  developing  many  particular  rate  differential 
equations.  It  should  be  noted  that  this  is  a  special  case  of  the  more  general  result. 

The  second  case  that  needs  to  be  developed  is  for  N  —  1.  As  previously  noted,  the 
(pseudo)  Saddle  Point  Expansion  is  valid  even  for  this  limiting  case,  so  we  may  return  to 
Equation  (85)  and  rewrite  it  as 


df(t)  =  A/ipi:1  (t)dt,  (91) 

=  A  fp(t)dt. 

We  note  immediately  that  the  jump  function  is  not  dependent  on  event  number  because  there 
is  only  one  event.  If  there  are  multiple  events,  however,  and  the  jump  process  is  repeated, 
this  differential  equation  becomes 


if  (t)  =  ^  A /<  n^p^(t)dt.  (92) 

n=  1 


Because  we  will  be  making  a  transformation  from  a  discrete  view  to  a  continuous  view,  the 
values  of  our  trajectories  will  be  transformed  from  integer  valued  to  non-integer  (nominally 
real)  valued.  It  is  thus  useful  to  examine  the  behavior  of  the  order  statistics  CDFs  when  the 
number  of  elements  in  the  collection  is  non-integer. 

B.  Extension  to  Non-Integer  Collections 

As  already  discussed,  Equation  (40),  the  probability  that  the  rth event  among  N  is  given 
by  a  simple  term  from  a  binomial  expansion, 

Pr:N  ( t )  =  ^  P  (I tY  [1  ~P  (t)]N-r,  (93) 

where  P  ( t )  is  the  probability  that  the  event  has  occurred  by  time  t  for  an  individual  element 
in  isolation.  These  probabilities  are  inherently  normalized, 

N  N  /  N\ 

E^m  =  E  L  Wru-p  (tr-r  p4) 

r= 0  r= 0  '  ' 

=  [P(t)  +  [l-P(t)]]N 

=  A 

from  their  basic  formulation,  so  this  expression  is  a  bit  of  a  tautology.  It  means,  however,  that 
so  long  as  N  has  integer  value,  the  collection  of  individual  event  probabilities  is  normalized 
for  all  values  of  time  regardless  of  any  evolution. 
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At  this  point,  we  lean  forward  and  anticipate  ourselves.  Specifically,  we  want  to  examine 
what  happens  when  we  allow  N  to  have  non-integer  (real)  values.  We  may  proceed  by 
considering  the  function 


ZR  =  [P(t)  +  \\-P  («)]]*;  (95) 

=  1, 

which  is  the  real  valued  equivalent  of  Equation  (38).  The  exponent  R  is  explicitly  non¬ 
integer.  For  convenience,  we  adopt  the  notation 

x  =  P  (t) ,  (96) 

V  =  [l--P(i)], 

so  that  we  may  write 

ZR  =  [x  +  y}R.  (97) 


This  may  appear  to  be  an  excessive  distinction,  but  it  may  be  recalled  that  this  is  exactly 
the  same  type  of  procedure  we  used  in  calculating  the  derivative  above.  It  is  also  a  common 
and  useful  distinction  in  manipulating  these  expressions  [Mohling  1982] . 


We  expand  this  function  in  a  Maclaurin  series  about  x  —  0,  giving  us  (in  a  naive  textbook 
fashion) , 


dkZ 


R 


k= 0 


dxk 


x=0 


X 

hi 


(98) 


The  derivatives  may  be  similarly  computed  as 


dkZR 

dxk 


x=0 


IK* 


U=o 


0 


y 


R-k 


(99) 


We  note  immediately  that  for  values  of  k  >  \  /?] ,  where  []  indicates  the  integer  value  of  the 
argument,  the  exponent  of  y  becomes  negative,  and  the  sign  of  the  product  begins  to  oscillate. 
The  series,  Equation  (98),  thus  becomes  oscillatory  beyond  this  value  of  k.  Further,  for  this 
part  of  the  series,  the  product  will  grow  in  magnitude  slower  than  A;!,  so  that  each  succeeding 
term  not  only  oscillates,  but  also  gets  smaller. 


This  series,  Equation  (98),  is  thus  Uniformly  Convergent  and  may  be  adequately  approx¬ 
imated  as 


[R] 


dkZ 


R 


k= 0 


dxk 


x=0 


X 

w 


(100) 


since  the  higher  terms  are  oscillatory  and  largely  cancel.  If  we  make  use  of  the  Gamma 
function,  we  may  rewrite  this  as 


ZR  ~ 


\ ^  r(A  + 1)  R-k^_ 

z-,oT(R+i^k)y  ki¬ 


rn) 
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This  is  the  real  number  equivalent  of  the  Binomial  Theorem.  In  a  like  manner  to  the  previous, 
we  may  write  the  PDF  of  the  kth  event  as 


Pk-.R 


T{R  +  1)  R_k  xk  1 

l  ift  -  1  -k)J  (fc  —  1) !’ 


(102) 


by  exactly  the  same  type  of  development  as  led  us  to  Equation  (41).  From  this,  we  may 
make  the  following  observations: 


1.  The  expected  value  of  the  probability  of  the  kth  of  R  events  is  well  defined  in  the  same 
manner  as  Equation  (43); 

2.  The  delta  function  behavior  of  the  PDF  of  the  kth  of  R  events  is  well  defined  in  the 
same  manner  as  Equations  (65  through  80);  and 

3.  The  summation  of  the  PDF  of  the  kth  of  R  events,  k  =  1..  [I?]  (actually  oo),  is  well 
defined  and  identical  in  form  to  Equations  (82  through  85).  This  latter  follows  from 
the  derivation  of  Stirling’s  Approximation  from  the  Gamma  function  with  the  factorial 
as  a  special  case. 


Thus,  this  approximation  not  only  has  components  that  look  like  the  purely  integer 
index  form  but  also  has  comparable  expansion  properties.  Based  on  these  observations,  the 
behavior  of  the  ”  Order  Statistics”  for  non-integer  and  integer  populations  is  sufficiently  close 
that  we  may  treat  them  identically.  The  results  derived  above  are  happily  applicable  to  real 
as  well  as  integer  numbers. 

C.  Deterministic  Rate  Equations 


It  may  be  recalled  that  meaningful  general  transformations  could  only  be  developed  for 
evenly  spaced  deterministic  event  times.  We  limit  our  immediate  discussion  to  such,  noting 
that  if  the  event  times  are  suitably  disordered,  or  have  some  pattern,  then  they  may  be 
represented  by  means  akin  to  those  just  developed  for  stochastic  events. 


For  single  sequences  (sets)  of  evenly  spaced  deterministic  event  times,  where  the  event 
jump  functions  are  independent  of  index,  we  identify  two  cases.  In  the  first  case,  events  begin 
in  index  zero,  that  is,  at  the  temporal  origin,  and  give  rise  to  a  rate  differential  equation, 


l  Am 

dt  At 


(103) 


where:  A f  (t)  =  the  jump  function,  and  1  /  At  =  the  continuous  representation  (approxima¬ 
tion)  of  the  jumps  that  occur  with  interval  At.  In  the  second  case,  the  jump  events  begin 
with  index  one,  and  give  rise  to  a  rate  differential  equation, 


l  Am 

dt/  At 


[1  _  e-2t/Ar  . 


(104) 


We  note  that  the  Equation  (104)  asymptotically  reduces  to  the  form  of  Equation  (103)  as 
time  increases. 
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D.  Comparison  of  Deterministic  and  Stochastic  Rate  Differential  Equations 

We  have  already  noted  that  for  stochastic  events,  the  rate  differential  equation  takes  the 
form 

^f(t)  =  &f(t)p(t)N(t),  (105) 

when  the  events  (jumps)  are  not  repetitive  and 

jtl  (t)  =  A/  (t)  p<R>  (t)  N  (t) ,  (106) 

when  the  events  are  repetitive.  In  some  cases,  deterministic  and  stochastic  rate  differential 
equations  may  take  the  same  form. 

The  first  case  of  this  is  the  most  obvious.  If  a  stochastic  process  is  uniformly  distributed 
on  [tl5 t2\,  then  pit)  =  1  / (t2  —  ti)  , ti  <  t  <  t2,  and  Equation  (105)  reduces  to 

jtf  ( t )  =  N(t),h<t<  t2,  (107) 

which  is  of  the  same  form  as  Equation  (103). 

Next,  if  a  repeated  process  (described  by  Renewal  Theory)  is  uniform  and  simple  with 
NED  distribution  p(t)  =  a  exp  (— at ) ,  a  —  1  /(t)  ,  (t)  —  /0°°  t'  p  (■ t ')  dt',  then  p ^  (t)  =  a , 
and  this  reduces  Equation  (106)  to 

=  (108) 

which  is  also  identical  in  form  to  Equation  (103).  Thus,  in  both  of  these  cases,  we  have 
stochastic  processes  with  rate  differential  equations  that  are  identical  to  deterministic  rate 
differential  equations. 

Lastly,  consider  a  repeated  process  that  is  alternating  in  two  NED,  p(t)  —  a  exp  (—at), 
and  q  (t)  =  (3  exp  (—/3t).  The  renewal  density  function  of  this  process  is 

p{R)  (t)  =  [1  -  e-(o+/3)4]  .  (109) 

This  allows  us  to  write  Equation  (106)  as 

jtf  (t)  =  A/  (t)N(t)[l-  e-<«+«‘]  ,  (110) 

which  is  of  the  same  form  as  Equation  (104). 

While  NED  renewals  are  essentially  the  simplest  that  may  be  formed,  they  are  nonetheless 
quite  common  because  of  the  frequency  of  occurrence  of  NED  processes  in  nature.  Regardless, 
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we  must  acknowledge  that  the  form  of  a  rate  differential  equation  does  not  uniquely  indicate 
that  the  process  it  represents  is  either  stochastic  or  deterministic. 

IX.  DETERMINISTIC  EXAMPLES 

We  now  essentially  change  our  approach  from  theory  to  example.  The  theoretical  basis 
of  Rate  Theory  having  been  established,  it  now  becomes  useful  to  consider  several  examples 
of  how  Rate  Theory  may  be  applied  to  the  development  of  Rate  Differential  Equations.  In 
this  first  section  of  examples,  we  concentrate  on  examples  of  deterministic  rate  processes, 
those  of  Type  1. 

A.  Pong 

Our  first  example  is  based  on  the  computer  game  Pong.  This  game,  one  of  the  earliest, 
is  actually  a  very  complicated  mechanics  problem.  In  essence,  it  is  the  Classical  (as  opposed 
to  Quantum)  Mechanical  one-dimensional  particle-in-the-box  problem. 

The  problem  consists  of  the  following:  a  one-dimensional  box  of  dimension  d  contains  a 
single  particle  moving  with  speed  of  magnitude  u0.  The  box  is  assumed  to  have  infinitely 
more  mass  than  the  particle.  Collisions  between  box  walls  are  elastic  and  lossless. 

As  a  result,  both  energy  and  momentum  are  conserved.  The  particle  moves  in  a  closed 
trajectory  (an  orbit  on  [0,  d])  between  the  two  walls  with  speed  (velocity)  of  ±Uo-  Nor¬ 
mally,  we  should  not  apply  Rate  Theory  to  mechanics  problems  since  the  formalisms  of 
Mechanics;  Newton’s  equations;  or  the  derivative  operations  of  the  Lagrangian,  Routhian, 
or  Hamiltonian  functions;  provide  a  direct  means  of  developing  the  differential  equations  of 
mechanics  [Goldstein  1950].  However,  because  this  problem  has  nonholonomic  constraints, 
it  is  useful  to  develop  a  Rate  Theory  example,  approximate  though  it  is,  to  show  a  special 
application  of  the  theory. 

Despite  the  nonholonomic  constraint,  this  problem  is  so  simple  that  an  algorithmic  so¬ 
lution  can  be  written  trivially.  If  we  define  the  transit  time  of  the  particle  from  one  side  of 
the  box  to  another  as  At,  then  we  have 


At  =  — .  (Ill) 

u0 

It  is  then  useful  to  define  the  integer  j*  as 

j*  =  integer  ,  (112) 

where  t  =  the  time  since  the  particle  satisfied  the  conditions  a;  (0)  =  0,u(0)  =  u0,  where 
x  (t) ,  u  (t)  are,  respectively,  the  position  and  speed  of  the  particle  at  time  t.  These  are  what 
we  might  consider  to  be  initial  conditions.  In  this  context,  j*  is  just  the  number  of  times 
that  the  particle  has  bounced  off  the  walls  since  the  initial  time. 
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From  this  we  may  write  that 


u0  (t  —  j*At) ,  j*  even 
d  —  Uq  (t  —  j*  At) ,  j*  odd  ’ 


(113) 


x(t)  = 

which  is  indeed  the  solution  of  the  motion  of  the  particle,  albeit  a  piecewise  solution. 
If  we  return  to  Equation  (44), 


i 

then  we  may  write  the  rate  differential  equation  of  the  particle’s  speed  as 

Jfu  =  J2Aui5  (t  -tj).  (114) 

3 

We  have  switched  indices  from  i  to  j  for  reasons  that  will  be  abundantly  clear  shortly. 

The  jump  times  tj  are  the  times  when  the  particle  collides  with  one  or  other  of  the  walls, 
which  is  just 

tj  —  jAt.  (115) 

When  the  collision  occurs,  the  speed  of  the  particle  switches  from  positive  to  negative,  so  we 
may  write  the  jump  function  in  a  still  general  form  as 

Auj  —  (— iy  Art,  (116) 

where  we  have  left  the  magnitude  of  the  jump  unspecified  to  allow  us  to  correct  for  the 
approximation.  This  is  a  case  where  the  jump  function  does  depend  on  index  but,  as  we 
shall  see,  is  simple  enough  (that  is,  linear)  that  we  may  treat  it  in  closed  form. 

Before  proceeding,  it  is  useful  to  recall  the  identity 

-1  =  e±ni,  (117) 

where  i  =  \/—l.  which  illuminates  our  need  to  change  index  notation.  This  identity  allows 
us  to  rewrite  Equation  (116)  as 

Auj  —  enijAu,  (118) 

which  in  turn  allows  us  to  write  the  jump  terms  as 

OO 

E  A'b 5  -  y)  =  E  e±nijAu8(t-jAt),  (119) 

3  3=1 

OO 

=  Au^e±nij8{t-jAt). 

3=1 
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We  may  now  take  the  Laplace  Transform  of  the  right  side  of  Equation  (119)  as 


A u  enij  S(t  —  jAt) 

3= i 

This  series  is  exactly  summable,  yielding 


3= 1 


e-sjAte±nij 


(120) 


-sjAt  ±nij 


3=1 


gsAt.^fni  _  ^  ’ 


(12!) 


which  we  now  expand  to  first  order  using  the  techniques  described  previously  for  evenly 
spaced  event  times  as 

1  , 

(122) 


£e 

3=1 


-sjAte±nij 

_ 


sAt  =1 =  7T  i 

At  this  point,  we  resolve  the  sign  ambiguity  by  using  each  one  equally 


£ 

3=1 


e~sjAte±irij 


1  1 

-  + 


2  \.sAt  +  7r i  sAt  —  nij 

sAt 


S2At 2  +  7T2' 


We  now  rewrite  this  as 


£e 

3=1 


-sjAte±nij 

r^j 


A^2  +  T 


(123) 


(124) 


which  we  immediately  (after  consulting  a  Laplace  Transform  Table)  recognize  as 


-l 


At  q2  + 

+  At2 


=  COS 


At  \AtJ 


Tit  \ 


- 


(125) 


This  allows  us  to  rewrite  Equation  (119)  as 


AU  ^  e±m{j+1)  S(t-  jAt)  ~  -  ^  cos  , 

3=1  '  V  '/ 


(126) 


and  thus  the  rate  differential  equation  becomes 


d  A  u  f  Tit  \ 

—u  ~  — —  cos  — 

dt  At  \At  J 


(127) 


Before  proceeding,  it  is  useful  to  replace  At  with  its  definition,  Equation  (111),  which 
now  gives  us 


d  A  u  Uq 

—u  ~ - - — 

dt  d 


cos 


(128) 
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We  now  apply  boundary  conditions  that  the  particle  start  at  one  side  of  the  box  initially 
(x  (0)  =  0)  and  be  at  the  opposite  side  of  the  box  an  appropriate  time  later  (x  (At)  =  d). 
This  gives  us 


Simple  elimination  gives  us 


and 


(132) 

(133) 


This  gives  us  a  speed  jump  factor  that  is  7r2  / 4  times  what  we  should  have  expected  based 
on  purely  mechanical  considerations.  The  resulting  approximate  rate  trajectory  is 

x(t)~^  l-cos(^)  ,  (134) 

which  we  compare  to  the  exact  trajectory,  Equation  (113),  in  Figure  9.  The  agreement  is 
quite  good,  largely  because  of  the  boundary  conditions  used. 

B.  Bouncing  Ball 

We  are  all  familiar  with  the  trajectory  of  a  ball  that  is  thrown  or  dropped,  follows  a 
ballistic  trajectory  impacts  the  ground,  bounces,  follows  another  ballistic  trajectory —  albeit 
usually  a  bit  smaller — impacts  and  bounces  again,  and  so  forth.  If  we  neglect  air  resistance, 
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Figure  9.  Comparison  of  Exact  and  Rate  Pong  Trajectories 


which  is  proportional  to  either  the  speed,  or  the  square  of  the  speed,  [Marion  1965]  and 
which  is  small  if  the  speed  of  the  ball  is  low,  then  the  trajectory  of  each  bounce  is  given  by 

x(t)  =  xQt,  (135) 

z  (t)  =  'o t  -  ^ gt 2, 

where:  x,  z  are  the  horizontal,  vertical  coordinates,  t  is  parametric  time,  *o  are  the  initial 
(t  —  0  measured  from  the  instant  when  the  bounce  commences)  velocity  components,  and  g 
is  the  acceleration  due  to  gravity.  The  problem,  of  course,  is  writing  component  differential 
equations  (vector  differential  equation)  for  the  entire  trajectory.  We  approach  this  problem 
using  our  rate  theory  formalism.  As  with  the  Pong  example,  both  jump  magnitudes  and 
times  are  deterministic. 

Most  discussion  of  collisions  in  elementary  textbooks  is  devoted  to  elastic  collisions  where 
energy  and  momentum  are  conserved.  As  we  know  from  observing  bouncing  ball  trajectories, 
the  collisions  of  the  ball  with  the  ground  are  inelastic  since  the  ball  eventually  comes  to  rest. 
Every  time  the  ball  hits  the  ground,  it  loses  energy.  (Unless  the  ball  is  made  out  of  Disney 
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Flubber® ,  in  which  case,  it  gains  energy;  regardless  the  formalism  is  the  same.)  Four 
mechanisms  may  be  identified  for  this: 

1.  The  ball  loses  a  fixed  increment  of  its  energy  with  each  collision; 

2.  The  ball  loses  a  fixed  fraction  of  its  energy  with  each  collision; 

3.  The  ball’s  momentum  is  decreased  by  a  fixed  increment  with  each  collision;  or 

4.  The  ball’s  momentum  is  decreased  by  a  fixed  fraction  with  each  collision. 


This  taxonomy,  while  simplistic,  follows  from  the  form  of  drag  loss.  We  now  consider 
each  of  these  losses  in  turn. 

1.  Fixed-Energy  Increment 

From  our  development  of  rate  theory  methodology,  we  may  write  velocity  component 
rate  equation  primitives  as 


XjS(t-Tj),  (136) 

3 

A  V  (*  “  ri)  “  9, 

3 

where:  A Xj,  A Zj  are  the  horizontal,  vertical  velocity  component  jumps  and  r:l  are  the  jump 
times.  The  jump  times  are  when  the  collisions  occur  and  are  obviously  the  same  for  both 
components. 

We  note  that  while  energy  and  momentum  are  not  conserved  during  the  collisions,  they 
are  conserved  during  the  bounces  since  we  are  neglecting  air  drag.  Further,  mass  may  be 
assumed  to  be  conserved.  If  we  adjust  the  origin  of  the  local  coordinate  system  to  each 
successive  collision  point,  we  may  place  the  potential  energy  zero  at  that  origin  so  long  as 
the  overall  collision  surface  is  level.  Accordingly,  we  may  write  the  energy  at  collision  as 


1 1  • 
—x 
dt 

d  • 

~rZ 

dt 


where  m  is  the  mass  of  the  ball. 


E 


m 

2~ 


(137) 


If  we  now  consider  this  in  a  circular  coordinate  system  (r,  6)  then  immediately  around 


the  collision, 


> 


rd 


,  so  we  may  approximately  reduce  Equation  (137)  as 


E  ~ 


m»  2 
— r  . 


2 


(138) 
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With  this  equation,  even  though  energy  and  thereby  momentum  are  not  conserved,  we 
may  take  momentum  as  proportionally  conserved.  That  is,  the  magnitude  of  momentum 
is  changed  in  the  collision,  but  the  ratios  of  the  momentum  components  are  taken  to  be 
conserved.  As  a  result,  if  we  designate  the  component  value  just  before,  after  the  collision 
with  subscripted  minus,  plus,  we  may  write  proportional  conservation  relations, 


x+ 

•  ~  •  > 

r_  r+ 

•  • 

•  • 

r_  r+ 

The  sign  change  reflects  the  ball  going  from  moving  down  to  up. 


We  may  invert  Equation  (138)  as 


(139) 


2  E 


m 


which  allows  us  to  rewrite  Equations  (139)  as 


x+  = 


=  / 


z+  =  —  / 


'E± 

Ej 

[eI 


E _ 


(140) 


(141) 


This  gives  us  a  means  of  calculating  the  jump  magnitudes. 

Before  proceeding  with  that,  however,  it  is  first  useful  to  number  the  bounces  and  des¬ 
ignate  the  energy  of  each  bounce  as  E} ,  j  =  0..,  where  the  zeroth  bounce  is  the  initial  one. 
By  Equation  (137),  the  initial  energy  is  just 

£0  =  -(xo  +zo  ),  (142) 

and  the  other  energies  are  just 

Ej  =  E0  —  jAE,  (143) 

where  A E  is  the  constant  energy  loss  per  bounce.  Quite  obviously,  the  ball  cannot  have 
negative  energy,  so  either  the  energy  loss  is  an  integer  multiple  of  the  initial  energy  or  the 
last  collision  reduces  the  energy  to  zero.  Either  case  defines  the  number  of  bounces. 


Since  we  are  neglecting  air  drag,  as  we  have  already  indicated,  the  ball’s  trajectory 
between  collisions  is  conservative.  From  this  we  know  that  the  horizontal  velocity  component 
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is  constant  over  the  trajectory  element,  and  the  vertical  velocity  component  at  the  beginning 
of  the  element  is  equal  in  magnitude  and  opposite  in  sign  of  the  component  at  the  end  of  the 
element.  Let  us  now  designate  the  velocity  components  at  the  beginning  of  each  bounce  or 
trajectory  element  as  xv  Zj.  With  this,  we  may  now  combine  the  momentum  proportionality 
equations,  Equations  (141)  with  these  definitions  of  energy  and  velocity  components  to  write 
equations  that  relate  each  bounce  to  the  preceding  one 


Xj  =  aq_i  < 


Z  j  Z  j — \  \ 


Ej-Z 

Ej- 1 


(144) 


Note  that  the  minus  sign  in  the  vertical  component  relationship  disappeared  because  of  the 
definition  of  the  components. 


If  we  iterate  these  relationships,  we  may  arrive  at  a  general  form 


Xj  =  x0 


Zj  =  Z  0 


En  ' 


(145) 


We  may  now  use  the  definition  of  the  element  energies,  Equation  (143)  to  rewrite  these  as 

(146) 


Xj  =  Xo 


•  I E0  —  jAE 


En 


Zj  =  Z  o 


.  E0-jAE 


En 


Since  we  are  ignoring  air  drag,  we  may  solve  for  the  collision  time  increments  using  the 
vertical  component  of  Equations  (135).  This  gives  us 


Atj 


2zj-i 

9 


We  rather  obviously  reduce  this  to 


At,- 


2zq 

9 


U~1)AE 

E0 


T 1 


U-i)ae 

TP  ’ 


(147) 


(148) 
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since  r0  =  0.  This  allows  us  to  calculate  the  jump  or  collision  times  as 


Tj 


J2Ark 

k=  1 


nE 


(k  -  1)  AE 


This  series  is  not  obviously  summable. 


(149) 


We  are  now  in  a  position  to  calculate  the  velocity  component  jumps.  The  horizontal 
component  jump  is  straightforward, 


A  *  *  * 

A  Xj  =  Xj  —  Xj- 1. 

This  can  be  slightly  elaborated  by  substituting  Equation  (146) 

A  #  * 

A  Xj  —  x0 

The  vertical  component  jump  is  similar, 

A  Zj  =  zj  +  Zj_  i, 


Eq  —  jAE  /  Eq  —  (j  —  1)  AE 


Eft 


Eft 


(150) 


(151) 


(152) 


except  for  the  sign  change  resulting  from  the  change  in  direction  of  vertical  motion  across 
the  collision.  Elaboration  follows  directly  and  we  shall  not  belabor  the  obvious. 


At  this  point,  we  close  our  discussion  of  this  example.  Fundamentally,  the  presence  of 
the  index  of  summation  under  the  radical  precludes  analytical  summation  of  the  series  in  the 
jump  differential  equations,  Equations  (136).  Since  analytical  forms  for  the  event  times  and 
jump  magnitudes  have  been  developed,  numerical  computation  of  these  is  straightforward, 
and  thereby  numerical  simulation  of  the  jump  differential  equations.  Elaboration  of  this  is 
outside  the  scope  of  this  presentation. 

2.  Fixed-Energy  Fraction 


We  now  change  to  a  loss  representation  where  the  ball’s  energy  is  decreased  by  a  fixed 
fraction  1  —  e  every  time  it  collides  with  the  ground.  This  permits  us  to  write  the  equivalent 
of  Equation  (143)  as 

Ej  =  e3  E{).  (153) 

From  this,  we  may  proceed  to  develop  the  jump  magnitudes  and  jump  times. 


Our  first  step  is  to  revise  the  relationships  for  the  initial  velocity  components  of  each 
bounce,  Equations  (145).  These  revisions  are 


•  •  j/2 

Xj  =  x0eJ/  , 


(154) 


Zj  =  z0e 


4/2 
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The  jump  magnitudes  then  follow  using  Equations  (150)  and  (152)  as 


and 


A*,-  =  zne^2 


(155) 


(156) 


The  jump  time  increments  then  follow  directly  from  Equations  (147)  and  (154)  as 

Atj  =  Tie^2 .  (157) 

This  gives  jump  times  via  Equation  (149)  as 


ri  =  ri  S  e 

k= 1 

which  we  may  re-index  as 

k= 0 

We  immediately  recognize  this  as  a  geometric  series,  permitting  us  to  sum  it  as 

1  -  eJ/2 

Tj  ~  Tl  1  —  e1/2  ’ 

giving  us  an  analytical  representation  of  the  jump  or  collision  times. 


:(fc-l)/2 


(158) 


(159) 


(160) 


Despite  this  advance  over  the  previous  case,  we  still  do  not  have  a  useful  rate  differential 
equation.  To  see  this,  we  take  the  jump  series  for  the  horizontal  component  in  Equation 
(136)  and  substitute  appropriately, 


(t  -  Ti)  =  (l  ~  4)  eJ/ 26  ~  Tl  |  _  gi/2  )  ’  (161) 

where  J*  is  the  number  of  the  last  collision  when  the  energy  goes  to  zero.  If  we  now  take 
the  Laplace  Transform  of  this  equation,  we  obtain 


£ 


-  j 


1  -  el/2  \ 

1  -  e1/2  )  ' 


(162) 


We  see  immediately  that  if  we  attempt  to  expand  the  exponential  in  the  manner  previously 
elaborated,  we  obtain  powers  of  \Z\112  in  the  denominator  of  the  expansion.  How  the  result¬ 
ing  series  may  be  summed  is  not  obvious.  Thus,  we  abandon  further  discussion  of  this  case, 
as  we  did  the  previous,  for  the  same  reason. 
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3.  Fixed-Momentum  Increment 


We  now  turn  our  attention  to  a  case  where  the  loss  of  energy  may  be  expressed  as  a  loss 
of  momentum  by  a  fixed  increment.  To  pursue  this,  we  return  to  Equation  (139)  and  make 
use  of  this  loss  of  momentum,  giving  us 


x+  =  X--^,  (163) 

r_ 

and  the  analogous  equation  for  the  vertical  component  of  velocity.  Using  our  previous 
definitions,  we  may  rewrite  this  equation  as 


and  iterating  as  before  produce 


Xj  =  Xj_i 


ri- 1 


Xj  =  x0- 


r0 


(164) 


(165) 


If,  following  the  manner  of  our  previous  definition,  we  define 

*j  =  r0  -  j  Ar, 

then  we  may  rewrite  Equation  (165)  as 


Xj  =  x0- 


.  r0  -  j  Ar 


r  o 


We  may  now  write  the  horizontal  jump  magnitudes  using  this  as 

* o  -  J  A*  r0  -  (j  -  1)  Ar 


Arc.  =  xq 


r  o 


r0 


Af 


=  -Xq- 


r  o 


This  is  pleasingly  simple.  The  vertical  jump  magnitudes  are  not  so  simple, 

Azj  =  z  o 


■r0  -  j  Ar  r0  -  (j  -  1)  Ar 


—  2^o 

which  at  least  is  linear  in  the  index. 


r  o  r  o 

ro  -  (.)  -  \)  Ar 


To 


(166) 


(167) 


(168) 


(169) 
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We  may  also  calculate  the  jump  time  increments  from  Equation  (147)  as 


A  Tj  =  T0 


=  T0 


r0  ~  (j  ~  1)  A r 
r0 

(j  ~  1)  A  r 


(170) 


To 


We  may  substitute  this  into  Equation  (149)  as 


T3  = 


Et° 

k= 1 


(k  —  1)  A r 


r  o 


(171) 


which  we  may  separate  and  partially  re-index  as 


Tj  =  T0 


3  *  •  j-1 

E^-E* 

,k=  1  ?0  k= 0  . 


(172) 


Both  of  these  series  are  summable  using  basic  finite  difference  methods,  giving  us 


Tj  =  T0 


.  _  A  r  j  ( j  -  1) 

J  9  C\ 

r0  Z 


(173) 


This  form  of  the  jump  or  collision  times  is  sufficiently  complex  that  it  prohibits  sum¬ 
mation  of  expansions  of  the  Laplace  Transform  of  the  jump  series  in  the  rate  differential 
equations.  Accordingly,  we  again  cease  at  efforts  short  of  simple  expressions  of  these  differ¬ 
ential  equations. 

4.  Fixed-Momentum  Fraction 

The  last  case  we  take  up  in  this  example  is  one  where  the  momentum  is  reduced  by  a 
fixed  fraction  1  —  /r  at  each  collision.  As  a  result,  we  may  form  the  initial  horizontal  velocity 
components  as 

Xj  =  xo/j ,  (174) 

with  corresponding  equation  for  the  vertical  components.  Since  we  have  performed  this  trail 
of  manipulations  three  times  previously,  we  have  omitted  several  obvious  steps. 

From  this,  the  horizontal  velocity  jumps  are 

A  *j  =  xqh3  ^1 - ^  ,  (175) 

while  the  vertical  velocity  jumps  are 


The  jump  time  increments  are 


A Tj  =  r0/T  \  (177) 

and  thence  the  jump  times  are 

~j  -  r°/,/,'~1-  (178) 

k= 1 

We  may  re-index  this  and  sum, 

j- 1 

ri  = 

k=0 

1  -  p? 

=  To- - , 

1  —  /i 

giving  us  a  form  essentially  the  same  as  we  obtained  for  the 
case. 

As  a  result,  we  again  have  a  situation  where  we  cannot  form  simple  rate  differential 
equations.  While  we  have  defined  the  problem  and  elaborated  a  solution  means,  reduction 
to  elegant,  compact  differential  equations  has  eluded  us  in  each  case.  Nonetheless,  we  may, 
in  each  case,  elaborate  the  primitive  rate  differential  equations  and  they  are  readily  amenable 
to  numeric  simulation. 

X.  STOCHASTIC  EXAMPLES 

We  now  consider  a  series  of  examples  where  the  jump  magnitude,  or  the  event  times,  or 
both  are  stochastic. 

A.  The  Encounter  Problem 

The  encounter  problem  is  one  that  is  common  to  many  disciplines  although  it  may  be 
called  by  another  name.  For  example,  in  Physics  of  War  I  referred  to  it  as  the  Target 
Acquisition  problem.  I  call  it  the  encounter  problem  here  because  several  disciplines  are 
concerned  with  meetings  between  two  entities  where  one  or  both  of  the  entities  are  moving, 
usually  randomly.  Alternately,  one  entity  may  be  moving  while  the  other  is  stationary,  and 
either  the  motion  is  random,  or  the  stationary  position  is  random,  or  both. 

Assume  that  we  are  concerned  with  some  bounded  region  of  size  A,  which  may  be  one-, 
two-,  or  three-dimensional.  Assume  that  A  is  broken  up  into  N  subregions  of  size  A.  Thus, 
NX  =  A.  Further  assume  that  A  is  sized  so  that  if  an  entity  is  present,  it  can  (but  not 
necessarily  must)  be  detected  (observed),  but  if  no  entity  is  present  and  no  detection  occurs 
in  a  reasonable  period  of  time,  progress  may  continue  with  reasonable  confidence.  For 
elaboration,  let  there  be  m  (<C  N )  entities  in  A,  and  assume  no  A  contains  more  than  one 
entity.  Further  assume  that  if  there  is  an  entity  present,  the  (discrete)  probability  it  will  be 
detected  is  pd- 


(179) 


constant  energy  loss  increment 


45 


We  may  treat  this  as  an  independent  search  (sampling)  of  the  N  subregions.  Thus,  the 
probability  that  an  entity  is  present  in  a  subregion  to  be  detected  on  any  sample  is  m  /N , 
and  the  independent  probability  that  an  entity  is  detected  on  a  sample  is  pcj  rri  / N .  From 
this,  we  may  write  the  probability  of  detecting  an  entity  on  the  first  sample  as 

777 

P(l)=Pdjjl  (180) 

the  probability  of  detecting  an  entity  on  the  second  sample  as 

p(2)^pd-[l-pd-\,  (181) 

which  is  just  the  probability  of  no  detection  by  the  first  sample  times  the  independent 
probability  of  detecting  an  entity  on  a  sample,  and  the  probability  of  detecting  an  entity  on 
the  third  sample  as 

P  (3)  =  Pd  J-j:  [l  Prf  J^\  ■  (182) 

From  this,  we  may  generalize  to  the  probability  of  detecting  an  entity  on  the  nth  sample  as 

777  r  777  ~\  n— 1 

P{n^PdN[l~PdN\  ■  (183) 

Now,  let  P  ( t )  be  the  continuous  probability  that  a  detection  has  not  occurred  by  time  t. 
Since  the  process  here  is  repeated  (and  presumed  stochastic),  the  archetype  rate  differential 
equation  is  given  by  Equation  (106), 

Pp  it)  =  (t)  v'RHt)  N  (t) , 

if  the  jump  function  is  independent  of  index.  In  this  case,  we  take  the  number  of  agents  to 
be  fixed  as  one,  reducing  this  equation  to 

±P(t)  =  AP(t)pW(t).  (184) 

If  we  now  examine  the  discrete  formulation  above,  Equations  (180  to  183),  we  see  that  we 
may  write  the  jump  function  as 

_  777  _ 

A  P(t)  =  -pd-P(t)r  (185) 

since  the  discrete  term  in  brackets  is  just  the  probability  an  entity  has  not  been  detected. 
As  a  result,  we  may  rewrite  the  rate  differential  equation  as 

jT5  (■ t )  =  ~Pd  ^ P  (■ t )  p(R)  (■ t ) .  (186) 
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We  note  immediately  that  if  the  repeated  process  is  NED  distributed  (or  is  deterministic 
with  equally  spaced  event  times,)  then  the  renewal  density  function  reduces  to  a  constant. 
Regardless,  we  may  solve  this  rate  differential  equation  trivially  by  inspection  as 

P(t)  =  1  -  exp  pd  ^  p{R)  (?)  dt'^j  .  (187) 

Since  the  quantity  we  are  commonly  interested  in  is  the  PDF,  we  obtain  that  by  differenti¬ 
ation  as 

Pencounter  (t)  =  Pd  (*)  eXp  Pd  (?)  dt'^j  ■  (188) 

As  before,  we  note  that  if  the  renewal  is  simple,  uniform  and  NED  distributed,  the  renewal 
density  function  is  a  constant.  If  the  single  sample  probability  of  detection,  the  number  of 
entities,  and  the  number  of  subregions  are  constant,  then  this  reduces  very  neatly  to 

Pencounter  (t)  =  Pd  P^  (*)  exP  (~Pd  ^  P ^  (*))  •  (189) 


If  we  define  the  quantity 


then  we  may  change  variables  and  rewrite  Equation  (188)  as 

Pencounter  (*^)  exp  (  x)  , 


from  which  we  may  immediately  calculate  that 

(R)  ,  \  _  i 

Pencounter  (%  )  !■ 


Reversing  the  change  of  variable,  we  obtain 

P{enlounter  (*)  =  Pd  ^  p(R)  (t)  , 


(190) 


(191) 


(192) 


(193) 


which  now  gives  us  all  the  information  we  need.  This  is  an  example  of  using  Rate  Theory 
to  convert  discrete  or  mixed  stochastic  processes  into  continuous  representation. 

B.  Drunkard’s  Walk  1 


A  dipsomaniac  is  standing  on  a  sidewalk.  Every  time  increment  At,  he  (she)  takes 
a  forward  step  of  length  Ax  with  probability  /,  or  a  backward  step  (same  length)  with 
probability  b,  or  stays  in  place  with  probability  s.  This  process  is  illustrated  in  Figure  10. 


By  inspection,  we  see  that  the  jump  does  not  depend  on  index  but  is  sufficiently  compli¬ 
cated  that  some  analysis  is  required.  The  appropriate  rate  differential  equation  is 


d_ 

Ft 


u  (x,  t ) 


A u  (x,  t ) 
At 


(194) 
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Figure  10.  Illustrative  One-Dimensional  Drunkard’s  Walk  Trajectories  and  Mean  Trajectory 


where:  u  (. x ,  t )  =  probability  the  dipsomaniac  is  at  position  x  at  time  t.  The  challenge  is 
now  to  form  the  jump  function  A u  (x,t).  Since  the  jump  is  presumed  to  be  instantaneous, 
we  may  make  use  of  the  same  continuous  representation  as  we  used  in  the  previous  example. 

Because  there  are  three  states  -  forward,  backward,  and  stay  in  place  -  we  have  to  consider 
four  quantities  for  any  position  x.  Prior  to  the  jump,  the  probability  the  dipsomaniac  is  at 
this  position  is  u(x,t).  There  is  also  probability  s  that  this  position  is  maintained.  There 
is  also  probability  /  that  the  dipsomaniac  will  move  from  x  to  x  +  Arc,  and  probability  b  of 
moving  to  x  —  Arc.  Recognizing  that  probability  flows  into  the  position  during  the  jump,  we 
may  write  the  before  jump  term  as 


^ before  U  (*£;  t)  , 


(195) 


and  the  after  term  as 


Uafter  =  /  U  (x  —  Ax,  t)  +  b  U  ( X  +  A.X,  t)  +  S  U  (rC,  t)  .  (196) 

Notice  that  the  forward  and  backward  probabilities  are  associated  with  positions  in  the 
opposite  directions.  This  follows  from  the  flow  of  probability  into  position  rc  from  position 
x  —  Ax  moving  forward  with  probability  /. 
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The  combination  of  these  two  give  us  the  jump  function  as 


Aw  ((£,  f  )  ^after  ^ before 

=  f  u(x  —  Ax,  t)  +  b  u(x  +  Ax,  t)  +  (s  —  1)  u  (x,  t ) 


(197) 


We  now  expand  the  right  side  of  this  equation  in  a  Taylor’s  Series  about  x.  This  gives  us 


Au(x,t)  =  f 


,  ,  d u(x,t)  .  d2u(x,t)Ax2  Tr^rr, 

u  (x,  t ) - Ax  +  A  ;  —  +  HOT 


dx 


dx 2 


+b 


,  ,  d  U(x,t)  .  d2u(x,t)Ax2 

u  (x,  t )  +  :  ’  ’Ax  +  ;  -  +  HOT 


dx 


dx 2 


+  (s  -  1)  u  (x,t) . 

From  conservation  of  probability,  we  may  simplify  this  somewhat  as 

a  i  \  /7  r.du(x,t).  d2u  (x,  t)  Ax2 

A  u(x,t)  =  {b-f)  ;  JAx  +  (b  +  f )  1  j 


dx 


dx2 


and  thence  the  rate  differential  equation 


d  .  .  .  A xdu(x,t)  (b  +  /)  Ax2  d2u  (x,  t) 

Wu(X,t)  =  (b-f)  ^ 


(198) 


(199) 


(200) 


We  see  immediately  that  this  is  the  differential  equation  associated  with  a  diffusion  process, 
the  first  right-side  term  being  associated  with  drift  and  the  second  with  diffusion  proper. 


C.  Drunkard’s  Walk  2 


The  previous  example  presumed  that  the  steps  were  taken  at  even,  deterministic  intervals. 
Now  we  elaborate  by  assuming  that  the  time  when  steps  occur  is  stochastic  with  PDF  p  (t) 
and  that  time  to  step  is  independent  of  step  length. 


In  this  case,  we  should  have  started  with  rate  differential  Equation  (106)  instead  of 
Equation  (103) 

^  u  (x,  t)  =  A u  (x,  t)  p(R)  (■ t ) .  (201) 

The  renewal  density  function  can  be  calculated  directly  since  we  know  the  PDF,  and  the 
jump  function  is  the  same  as  developed  immediately  previous.  As  a  result,  we  may  write  the 
rate  differential  equation  as 


d_ 

Ft 


u  (x,  t ) 


(6-/)  A 

+AHl^pw(t)d2u^ 


dx2 


(202) 


which  we  still  recognize  as  a  diffusion  differential  equation.  As  before,  if  the  time  to  step 
PDF  is  NED,  then  this  rate  differential  equation  reduces  back  to  the  form  of  Equation  (200). 
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D.  Drunkard’s  Walk  3 


We  now  elaborate  when  the  step  length  is  still  independent  of  time  to  step  but  becomes 
stochastic  with  PDF  q  (x).  For  simplicity,  we  assume  that  forward  and  backward  steps  are 
identically  distributed.  Elaboration  to  either  different  distributions  or  a  single  distribution 
for  forward  and  backward  steps  is  straightforward  and  we  shall  not  treat  it  here. 


Our  starting  point  is  Equation  (201),  and  we  proceed  through  the  same  process  of  devel¬ 
oping  the  jump  function  as  before.  This  gives  us  a  rate  differential  equation 


d_ 

Ft 


U  ( X ,  t) 


(o SM;M) 


dx 2 


which  is  the  same  as  Equation  (202)  except  that  we  now  interpret  Ax  as  a  random  variable. 
Accordingly,  we  need  to  calculate  the  expected  value  of  this  equation  with  respect  to  this 
random  variable.  This  is  straightforward  and  yields 

Jh  (x,  t)  =  (b  -  /)  (Ax)  plR>  ( t )  9UgXJ  —  (203) 

+h±h(A^>  ^>  (*>%*>, 

where  (Arc) ,  (Ax2)  are  the  expected  values  of  the  step  length,  step  length-squared.  The 
left  side  is  not  changed  since  it  is  independent  of  step  size.  The  definition  of  the  standard 
deviation  allows  us  to  rewrite  this  as 


d_ 

Ft 


u  (x,  t) 


(!»  -  /)  (Ax) 

+  h±/)(ff2  +  (Ax)2)  p™(t) 


d2u  (rc,  t) 
dx 2 


which  illustrates  a  difference  from  what  is  normally  presented  in  textbooks. 


(204) 


E.  The  Logistics  Equation 


While  this  section  and  its  differential  equation  are  called  the  Logistics  Equation,  I  am 
going  to  approach  it  from  the  standpoint  of  the  spread  of  disease  because  I  find  it  easier  to 
understand  and  develop  in  that  context. 

The  disease  in  question  is  non-lethal  -  so  we  don’t  have  to  worry  about  counting  those 
who  die  -  and  is  spread  only  through  direct  contact,  which  may  be  either  physical  contact  or 
near  proximity.  The  quantity  d  (t)  is  the  number  of  people  at  time  t  who  have  the  disease. 
We  take  the  total  population,  which  is  fixed,  to  be  N  and  the  likelihood  that  the  disease  is 
caught,  given  contact,  is  l.  Obviously,  we  also  assume  that  no  one  ever  recovers. 

All  of  these  restrictions  can  be  relaxed  or  eliminated,  but  we  are  aiming  toward  a  result 
here  and  don’t  want  to  get  too  complicated  too  fast. 
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The  starting  point  is  Equation  (106), 

^-d  (t)  =  Ad  (t)  p^  ( t )  N  (t) .  (205) 

LLL 

We  note  immediately  that  the  number  of  agents  causing  events  is  just  the  number  of  people 
with  the  disease,  so  we  may  rewrite  this  immediately  as 

^-d  (t)  =  Ad  ( t )  p ^  ( t )  d  (t) .  (206) 

If  an  encounter  (event)  occurs  between  two  diseased  people,  the  jump  is  zero,  since  both 
have  the  disease;  if  an  encounter  occurs  between  a  diseased  person  and  an  undiseased  person, 
then  the  undiseased  person  may  contract  the  disease  with  likelihood  l.  Thus,  we  need  only 
count  encounters  between  diseased  and  undiseased  people  and  the  jump  is  stochastic  with 
magnitude  l. 

The  renewal  density  function  represent  the  encounters,  and  since  we  have  a  zero  jump 
for  encounters  between  two  diseased  people,  we  only  have  to  count  encounters  between  a 
diseased  and  an  undiseased  individual.  This  is  a  variant  on  the  Encounter  Problem,  so  we 
may  use  Equation  (193), 

pinlounter  (*)  =  Pd  ^  q{R)  (t)  , 

to  build  the  renewal  density  function  for  the  Logistics  Equation.  In  this  case,  the  number 
of  subregions  is  just  the  total  number  of  people,  which  is  N  —>■  N,  the  number  of  entities 
that  can  be  encountered  is  the  number  of  undiseased  people,  which  is  m  — >  N  —  d  (t) ,  the 
probability  of  detection  is  one  -  pd  — > ►  l,4  and  the  sample  renewal  density  rate  is  just  the  rate 
of  samples,  so  r/rd  (t)  — >  (t). 

As  a  result,  the  rate  differential  equation  becomes 

|d(t)  =  iiV~Af(*)  qW(t)d(t).  (207) 

This  equation  may  be  rearranged  as 

A  ]  n{R)  (f\ 

—d  (t)  =  l  q(R)  (' t )  d  (■ t )  -  N  d2  (■ t ) ,  (208) 

which  is  a  form  commonly  expressed  as  the  Logistics  Equation  [Ball  1985]. 

As  a  matter  of  clarification,  in  this  model,  the  initial  condition  has  to  be  d  (to)  >  0. 
Otherwise,  no  one  ever  gets  the  disease. 

4  Equivalently,  we  could  have  counted  the  encounters  not  as  such  but  as  disease  transmission.  In  that 
case,  the  probability  of  detection  would  be  the  likelihood  of  contracting  the  disease  and  the  jump  function 
one. 
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(209) 


F.  Predator  and  Prey 


The  Lotka-Volterra  equations,  [Ball  1985] 


dx 


—  =  Ax  —  Bxy, 
dt 

^  =  —  Cx  +  Dxy, 
dt 


are  classical  differential  equations  describing  the  dynamics  of  a  population  x  (t)  of  prey  and 
a  population  y  (t)  of  predators.  The  first  right-side  term  in  the  prey  equation  represents  the 
increase  in  the  prey  population,  and  the  second  term  its  decrease.  The  forms  are  reversed 
in  the  predator  equation. 


We  assume  that  these  prey  and  predator  constitute  a  closed  system  and  that  there  is 
sufficient  food  for  the  prey  that  they  never  lack. 


In  a  closed  system  where  members  do  not  migrate,  an  increase  in  the  number  of  members 
is  the  result  of  reproduction.  Reproduction  may  take  two  forms  (that  we  know  of):  sexual 
and  asexual.  In  asexual  reproduction,  one  parent  divides  into  two  or  more  offspring.  In 
sexual  reproduction,  two  parents  must  mate  to  produce  one  or  more  offspring,  n. 

In  terms  of  rates,  asexual  reproduction  produces  n  offspring,  where  n  may  be  the  ex¬ 
pected  value  of  some  random  variable,  so  we  have  a  rate  of  increase  of  asexually  reproducing 
population  as 

77/  —  1 

r asexual  ~  ~^TN  (0  ’  (210) 

where:  At  =  the  time  between  reproductions  (assuming  these  times  are  deterministic  or 
stochastic  with  NED  distribution,  otherwise,  =  1  / (t) ,  s ^  ( t )  =  renewal  density  function 
of  reproductions),  and  N  (t)  =  (reproducible)  population  at  time  t.  If,  on  the  other  hand, 
the  reproduction  is  sexual,  then  the  rate  of  that  reproduction  should  be 

77/  o 

5e»I~^(t) ,  (211) 

on  the  assumption  that  the  population  is  evenly  divided  between  male  and  female  members. 
(Whence  the  factor  of  2,  since  only  half  of  possible  encounters  will  be  between  opposite  sex 
individuals.)  The  form  of  this  rate  follows  directly  from  the  encounter  problem. 

There  are  two  major  causes  of  death  for  prey  and  one  for  predators.  Both  may  die  of  old 
age,  which  has  a  rate  of 

r  death - ^V(i),  (212) 

where:  At  =  the  average  life  span  (again  either  deterministic  or  NED).  This  assumes  we  have 
not  been  keeping  up  with  ages  and  death  by  old  age  occurs  uniformly  among  the  population 
(which  is  a  bit  of  a  misnomer.)  Death  by  being  eaten,  which  applies  uniquely  to  prey,  has 
the  rate 

r eaten - N  (*)  M  (*)  »  (213) 
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where:  At  =  the  encounter  time  for  predator  and  prey,  p  =  the  probability  that  a  predator 
will  catch  and  eat  a  prey  (actually  kill  a  prey)  in  an  encounter,  and  N  (t) ,  M  (t)  =  the 
populations  of  prey  and  predator. 


If  we  now  assume  that  both  predator  and  prey  reproduce  asexually,  then  the  predator- 
prey  rate  differential  equations  should  be 


dx 

dt 

dy 

dt 


n  —  1 


At. 

n 


-x  — 


x  — 


P 


repro 

'  1 


At: 


/ 

repro 


-y 


At 0id  age 

1 

At' 

old  age 


At 


-xy, 


eaten 


(214) 


since  eating  prey  does  not  cause  new  predators  directly  to  be  produced.  Note  that  the 
predator  quantities  are  primed  here.  If  both  predators  and  prey  reproduce  sexually,  then 
the  rate  differential  equations  may  be 


dx 

dt 

dy 

dt 


n 


2A  t. 
n' 


-x2  — 


repro 

/ 


At 0id  age 

1 


X  — 


V 


At 


-xy, 


eaten 


2At 


i 

repro 


-r 


At 


-y- 


old  age 


(215) 


Notice  that  neither  set  of  these  rate  differential  equations  look  like  the  Lotka-Volterra  equa¬ 
tions. 


We  might  want  to  account  for  a  predator  starving  to  death  if  it  doesn’t  get  to  eat  a  prey 
in  a  certain  amount  of  time.  If  we  say  that  a  predator  will  starve  if  it  goes  k  encounters  with 
a  prey  without  eating.  If  the  time  between  encounters  is  A tenCounter  and  the  probability  of 
eating  a  prey  on  an  encounter  is  ~  peatx  /£,  where  £  =  the  ratio  of  total  hunting  area  to 
encounter  area,  then  the  rate  of  starvation  is  approximately 


T* starve  ^ 


1 


kAt 


encounter 


PeatX 

£ 


y, 


(216) 


and  the  sexual  rate  differential  equations  become 


dx 

dt 

dy 

dt 


n 


2A  t 

n 


-x 


repro 

/ 


2A  t 
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repro 


l 


-x 
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At  old  age  At 

1 

—y 


-xy, 


eaten 

1 


At 
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kAt, 


encounter 


PeatX 

£ 
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(217) 


If  peatx  /£  <  1,  then  we  have  the  forms 


dx 

dt 

dy 

dt 


n 


2A  t 


-x2  — 


n 


repro 

/ 


2At 


/ 

repro 


- y 2- 


Ataid  age 

1 


X  — 


P 


At  ~XVl 

eaten 

-1  _  kpeatX 

-L 


At 


-y  - 


old  age 


kAt, 


y f 


encounter 


(218) 
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which  gets  us  closer  to  the  form  of  the  Lotka-Volterra  equations,  especially  for  the  asexual 
reproduction  forms. 

This  is  a  rather  obvious  example  where  rate  theory  produces  something  strikingly  different 
from  what  is  commonly  used  on  a  more  ad  hoc  basis. 

G.  Lanchester  Equations 

The  Lanchester  Equations  [Lanchester  1916]  describe  losses  to  a  force  in  combat  due  to 
the  exchange  of  fire.  An  elaboration  of  these  equations  is  commonly  used  to  describe  forces 
comprised  of  more  than  one  type  of  combat  element,  but  we  shall  restrict  our  discussion  here 
to  the  homogeneous  problem.  Considerable  discussion  of  the  validity  of  these  differential 
equations  rages,  in  many  cases  because  of  the  lack  of  understanding  of  what  they  represent. 
We  derive  these  equations  to  their  commonly  used  (and  restricted)  forms. 

The  Lanchester  Equations  are  normally  paired  and  symmetric,  although  this  is  not  nec¬ 
essary  and  is  more  convenience  and  tradition  than  necessity.  Additionally,  while  they  admit 
(at  least  for  one  form)  changes  in  “force  strength"  due  to  reinforcement,  they  are  commonly 
started  impulsively  while  real  combat  is  seldom  so  ordered  [Fowler  2004], 

We  consider  two  cases  where  each  individual  weapons’  fire  has  a  lethality  area/volume 
that  is  less  than  that  of  a  single  combat  element  or  more  than  that  of  one.  Thus,  we 
distinguish  between  situations  of  point  lethality  where  at  most  one  combat  element  may  be 
attrited  and  area  (volume)  lethality  where  more  than  one  combat  element  may  be  attrited. 

We  label  the  two  sides  as  enemy  and  friendly,  and  designate  the  number  of  elements 
at  time  t  as  E  (t)  for  the  enemy  force,  and  F  (t)  for  the  friendly  force.  We  also  limit  our 
discussion  to  only  one  side’s  suffered  attrition. 

1.  Point  Lethality 

Point  lethality  is  commonly  associated  with  direct,  that  is,  line-of-sight,  fire.  It  commonly 
consists  of  two  processes:  finding  an  enemy  element  to  shoot  at  and  the  actual  shooting 
itself.  In  most  cases,  the  first  is  a  variant  of  the  encounter  problem.  We  can  approximate 
the  target-finding  process  with  a  PDF  ala  Equation  (188) 

P encounter  (t)  =  Pd  ^  P^  (t)  exp  J  pd  ^  p(R)  (t')  dt'^j  (219) 

E(t)  (  E{t)\ 

-  exp(^-a—  tj, 

where:  A  (t)  =  the  area  (or  volume,  depending  on  the  search  and  sensor  forms,)  occupied  by 
the  enemy  force,  the  search  process  is  assumed  to  be  uniform,  simple,  and  NED,  and  the  ratio 
E  (t)  /A  (t)  is  assumed  to  be  constant  or  slowly  varying.  Since  we  will  be  leading  up  to  some 
limiting  behaviors,  the  latter  is  not  critical  to  the  development  of  Lanchester’s  Equations 
but  is  in  the  general  context  of  attrition  modeling.  Also,  because  of  these  restrictions  and 


54 


assumptions,  all  other  factors  may  be  lumped  in  the  coefficient  a.  Note  that  we  may  develop 
A  (t)  from  the  definition  of  A  and  A  in  our  discussion  of  the  encounter  problem  above. 


The  shooting  process  is  assumed  to  be  NED  with  rate  (3.  As  a  result,  we  may  write  the 
point-lethality  process  renewal  density  function  as 


Plaint  (*) 


a 


m  a 

A(t)P 


01  W)  +  P 


E(t) 

m 


-p)t 


(220) 


from  Equation  (109).  (We  are  considering  here  the  case  of  friendly  elements  shooting  at 
enemy  elements.) 


Since  this  is  point  lethality,  the  jump  is  uniquely  —1,  since  only  one  element  is  attrited 
and  the  attrition  reduces  the  force  strength.  Thus,  from  Equation  (106)  we  may  write  the 
rate  differential  equation  of  the  enemy  force  strength  as 


d 

dt 


E(t) 


a 


m  a 

A(t)P 


a 


Mil  _i_  r  l 

A(t)  ^  P 


E(t) 

Ait) 


F{t). 


(221) 


We  now  wish  away  the  time  dependent  term  in  the  brackets.  Two  related  but  different 
arguments  permit  this.  The  first,  which  is  similar  to  that  used  in  Bonder-Farrell  Theory 
[Fowler  2004,  and  references  therein]  is  that  the  rate  is  defined  ala  Blackwell’s  Theorem  in 
the  infinite  time  limit.  The  second  is  that  we  are  interested  in  the  steady  state  or  equilibrium 
behavior.  Regardless,  Equation  (221)  simplifies  to 


d 

dt 


E(t ) 


a 


m  a 

A(t)P 
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M0  I  r 
A(t)  ^  P 


F{t). 


(222) 


We  now  assume  that  the  time  to  find  a  target  is  much  less  than  the  time  to  shoot  and 
attrit  a  target.  This  is  just  an  investigation  of  limiting  behavior.  As  a  result, 

<*!§»  A  (223) 

and  Equation  (222)  may  be  reduced  to 

jE(t)  =  -PF(t),  (224) 

which  is  a  linear  differential  equation,  but  in  Lanchester  Theory  terms  it  is  called  a  quadratic 
Lanchester  equation  (because  of  the  form  of  it  and  its  symmetric  complement’s  normal 
form  solution).  This  demonstrates  that  point  lethality  reduces  to  a  quadratic  Lanchester 
differential  equation  when  target  density  is  high,  target  shooting  is  slow,  or  target  finding  is 
fast. 
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When  the  opposite  is  true,  and  the  time  to  find  a  target  is  much  greater  than  the  time 
to  shoot  and  attrit  a  target,  Equation  (222)  reduces  to 

(226) 

Now,  if  the  density  of  enemy  elements  stays  fixed  over  time,  this  further  reduces  to 

5BW  =  -“fSF(th  (226) 

which  is  again  a  quadratic  Lanchester  differential  equation.  If  the  target  density  varies,  as 
would  be  the  case  if  the  area  occupied  stays  fixed,  then  Equation  (222)  reduces  to 

iE^=-amF(t)’  <227) 

which  is  a  non-linear  rate  differential  equation,  but  in  Lanchester  terms  it  is  a  linear  Lanchester 
differential  equation.  This  demonstrates  the  other  extreme  behavior. 

2.  Area  Lethality 

Area  lethality  fire  is  normally  directed  since  the  weapons  do  not  have  line-of-sight  to 
their  targets.  In  such  circumstances  and  assuming  that  the  firing  is  simple,  uniform,  and 
NED  with  rate  7,  we  may  again  use  Equation  (106) 

(t)  —  AE/y  F  (t) .  (228) 

The  jump  is  the  number  of  enemy  elements  killed  per  friendly  shot.  This  should  be 

(approximately) 

E  (t) 

AE  ~  —  Pk\s  A  7X  a  Lethal,  (229) 

where:  pk\s  =  the  probability  of  kill  in  the  lethal  area  given  a  shot,  a  Lethal  =  lethal  area  per 
shot,  and  E  (t)  /A  (t)  =  the  density  of  enemy  elements.  Substitution  of  this  yields 

jtE  (t)  =  ~Pk\s  jjlaLethan  F  (t)  ■  (230) 


If  we  now  take  the  limit  that  the  density  stays  constant,  then  this  reduces  to 

JtE  0)  =  -Pk\s  ^EaLethall  F  it) ,  (231) 

which  is  a  quadratic  Lanchester  differential  equation.  If  the  occupied  area  stays  fixed,  then 
the  limit  is 

jtE  (*)  =  ~Pk\s  ^ EaLethall  E  (t) ,  (232) 
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which  is  a  linear  Lanchester  differential  equation. 

Regardless  of  whether  we  are  considering  point  or  area  lethality  in  the  long  time  or  equi¬ 
librium  limit,  we  may  obtain  either  a  linear  or  a  quadratic  Lanchester  differential  equation, 
depending  on  limiting  assumptions. 
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